[u/mrichter/AliRoot.git] / PWGGA / EMCALJetTasks / AliJetResponseMaker.cxx

// $Id: AliJetResponseMaker.cxx  $
//
// 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() : 
  AliAnalysisTaskEmcal("AliJetResponseMaker"),
  fMCPartName("Tracks"),
  fMCJetsName("mcJets"),
  fMCParts(0),
  fMCJets(0),
  fHistMCJetPhiEta(0),
  fHistMCJetsPt(0),
  fHistMCJetsPtTrack(0),
  fHistMCJetsPtClus(0),
  fHistMCJetsPtNonBias(0),
  fHistMCLeadingJetPt(0),
  fHistMCJetsNEFvsPt(0),
  fHistMCJetsZvsPt(0),
  fHistJetPhiEta(0),
  fHistJetsPt(0),
  fHistJetsPtTrack(0),
  fHistJetsPtClus(0),
  fHistJetsPtNonBias(0),
  fHistLeadingJetPt(0),
  fHistJetsNEFvsPt(0),
  fHistJetsZvsPt(0)
{
  // Default constructor.
}

//________________________________________________________________________
AliJetResponseMaker::AliJetResponseMaker(const char *name) : 
  AliAnalysisTaskEmcal(name),
  fMCPartName("Tracks"),
  fMCJetsName("mcJets"),
  fMCParts(0),
  fMCJets(0),
  fHistMCJetPhiEta(0),
  fHistMCJetsPt(0),
  fHistMCJetsPtTrack(0),
  fHistMCJetsPtClus(0),
  fHistMCJetsPtNonBias(0),
  fHistMCLeadingJetPt(0),
  fHistMCJetsNEFvsPt(0),
  fHistMCJetsZvsPt(0),
  fHistJetPhiEta(0),
  fHistJetsPt(0),
  fHistJetsPtTrack(0),
  fHistJetsPtClus(0),
  fHistJetsPtNonBias(0),
  fHistLeadingJetPt(0),
  fHistJetsNEFvsPt(0),
  fHistJetsZvsPt(0)
{
  // Standard constructor.

  fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
}

//________________________________________________________________________
AliJetResponseMaker::~AliJetResponseMaker()
{
  // Destructor
}

//________________________________________________________________________
void AliJetResponseMaker::UserCreateOutputObjects()
{
  // Create user objects.

  OpenFile(1);
  fOutput = new TList();
  fOutput->SetOwner();

  fHistJetsPt = new TH1F("fHistJetsPt", "fHistJetsPt", fNbins, fMinPt, fMaxPt);
  fHistJetsPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
  fHistJetsPt->GetYaxis()->SetTitle("counts");
  fOutput->Add(fHistJetsPt);
  
  if (fAnaType == kEMCAL) {
    fHistJetsPtClus = new TH1F("fHistJetsPtClus", "fHistJetsPtClus", fNbins, fMinPt, fMaxPt);
    fHistJetsPtClus->GetXaxis()->SetTitle("p_{T} [GeV/c]");
    fHistJetsPtClus->GetYaxis()->SetTitle("counts");
    fOutput->Add(fHistJetsPtClus);
  }
  
  fHistJetsPtTrack = new TH1F("fHistJetsPtTrack", "fHistJetsPtTrack", fNbins, fMinPt, fMaxPt);
  fHistJetsPtTrack->GetXaxis()->SetTitle("p_{T} [GeV/c]");
  fHistJetsPtTrack->GetYaxis()->SetTitle("counts");
  fOutput->Add(fHistJetsPtTrack);
  
  fHistJetsPtNonBias = new TH1F("fHistJetsPtNonBias", "fHistJetsPtNonBias", fNbins, fMinPt, fMaxPt);
  fHistJetsPtNonBias->GetXaxis()->SetTitle("p_{T} [GeV/c]");
  fHistJetsPtNonBias->GetYaxis()->SetTitle("counts");
  fOutput->Add(fHistJetsPtNonBias);
  
  fHistLeadingJetPt = new TH1F("fHistLeadingJetPt", "fHistLeadingJetPt", fNbins, fMinPt, fMaxPt);
  fHistLeadingJetPt->GetXaxis()->SetTitle("p_{T} [GeV]");
  fHistLeadingJetPt->GetYaxis()->SetTitle("counts");
  fOutput->Add(fHistLeadingJetPt);
  
  fHistJetsZvsPt = new TH2F("fHistJetsZvsPt", "fHistJetsZvsPt", fNbins, 0, 1.2, fNbins, fMinPt, fMaxPt);
  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, fMinPt, fMaxPt);
    fHistJetsNEFvsPt->GetXaxis()->SetTitle("NEF");
    fHistJetsNEFvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
    fOutput->Add(fHistJetsNEFvsPt);
  }
}

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

  Int_t maxJetIndex  = -1;
  Int_t maxMCJetIndex  = -1;

  DoJetLoop(maxJetIndex, fJets, fTracks, fCaloClusters); 
  if (maxJetIndex < 0) 
    return;
  AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(fJets->At(maxJetIndex));
  if (!jet) 
    return;

  DoJetLoop(maxMCJetIndex, fMCJets, fMCParts);  
  if (maxMCJetIndex < 0) 
    return;
  AliEmcalJet* mcjet = dynamic_cast<AliEmcalJet*>(fMCJets->At(maxMCJetIndex));
  if (!mcjet) 
    return;

  fHistLeadingJetPt->Fill(jet->Pt());
  fHistMCLeadingJetPt->Fill(mcjet->Pt());
}

//________________________________________________________________________
void AliJetResponseMaker::DoJetLoop(Int_t &maxJetIndex, TClonesArray *jets, TClonesArray *tracks, TClonesArray *clusters)
{
  // Do the jet loop.

  if (!fJets)
    return;

  Int_t njets = jets->GetEntriesFast();

  Float_t maxJetPt = 0;
  for (Int_t ij = 0; ij < njets; ij++) {

    AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(jets->At(ij));

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

    if (!AcceptJet(jet))
      continue;

    fHistJetsPtNonBias->Fill(jet->Pt());

    if (jet->MaxTrackPt() > fPtBiasJetTrack)
      fHistJetsPtTrack->Fill(jet->Pt());
    
    if (fAnaType == kEMCAL && jet->MaxClusterPt() > fPtBiasJetClus) 
      fHistJetsPtClus->Fill(jet->Pt());
    
    if (jet->MaxTrackPt() < fPtBiasJetTrack && (fAnaType == kTPC || jet->MaxClusterPt() < fPtBiasJetClus))
	continue;

    if (maxJetPt < jet->Pt()) {
      maxJetPt = jet->Pt();
      maxJetIndex = ij;
    }

    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, tracks);
      if (track)
	fHistJetsZvsPt->Fill(track->Pt() / jet->Pt(), jet->Pt());
    }

    if (fAnaType != kEMCAL || !clusters)
      continue;

    for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
      AliVCluster *cluster = jet->ClusterAt(ic, clusters);

      if (cluster) {
	TLorentzVector nPart;
	cluster->GetMomentum(nPart, fVertex);
	fHistJetsZvsPt->Fill(nPart.Et() / jet->Pt(), jet->Pt());
      }
    }
  }
}

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

  AliAnalysisTaskEmcal::RetrieveEventObjects();
  
  if (strcmp(fMCJetsName,"")) {
    fMCJets = dynamic_cast<TClonesArray*>(MCEvent()->FindListObject(fMCJetsName));
    if (!fMCJets) {
      AliWarning(Form("Could not retrieve MC jets %s!", fMCJetsName.Data())); 
    }
  }
}

//________________________________________________________________________
void AliJetResponseMaker::Terminate(Option_t *) 
{
  // Called once at the end of the analysis.
}
Commit	Line	Data
2949a2ec	1	// $Id: AliJetResponseMaker.cxx $
	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() :
	28	AliAnalysisTaskEmcal("AliJetResponseMaker"),
	29	fMCPartName("Tracks"),
	30	fMCJetsName("mcJets"),
	31	fMCParts(0),
	32	fMCJets(0),
	33	fHistMCJetPhiEta(0),
	34	fHistMCJetsPt(0),
	35	fHistMCJetsPtTrack(0),
	36	fHistMCJetsPtClus(0),
	37	fHistMCJetsPtNonBias(0),
	38	fHistMCLeadingJetPt(0),
	39	fHistMCJetsNEFvsPt(0),
	40	fHistMCJetsZvsPt(0),
	41	fHistJetPhiEta(0),
	42	fHistJetsPt(0),
	43	fHistJetsPtTrack(0),
	44	fHistJetsPtClus(0),
	45	fHistJetsPtNonBias(0),
	46	fHistLeadingJetPt(0),
	47	fHistJetsNEFvsPt(0),
	48	fHistJetsZvsPt(0)
	49	{
	50	// Default constructor.
	51	}
	52
	53	//________________________________________________________________________
	54	AliJetResponseMaker::AliJetResponseMaker(const char *name) :
	55	AliAnalysisTaskEmcal(name),
	56	fMCPartName("Tracks"),
	57	fMCJetsName("mcJets"),
	58	fMCParts(0),
	59	fMCJets(0),
	60	fHistMCJetPhiEta(0),
	61	fHistMCJetsPt(0),
	62	fHistMCJetsPtTrack(0),
	63	fHistMCJetsPtClus(0),
	64	fHistMCJetsPtNonBias(0),
65	fHistMCLeadingJetPt(0),
66	fHistMCJetsNEFvsPt(0),
67	fHistMCJetsZvsPt(0),
68	fHistJetPhiEta(0),
69	fHistJetsPt(0),
70	fHistJetsPtTrack(0),
71	fHistJetsPtClus(0),
72	fHistJetsPtNonBias(0),
73	fHistLeadingJetPt(0),
74	fHistJetsNEFvsPt(0),
75	fHistJetsZvsPt(0)
76	{
77	// Standard constructor.
78
79	fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
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
97	fHistJetsPt = new TH1F("fHistJetsPt", "fHistJetsPt", fNbins, fMinPt, fMaxPt);
98	fHistJetsPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
99	fHistJetsPt->GetYaxis()->SetTitle("counts");
100	fOutput->Add(fHistJetsPt);
101
102	if (fAnaType == kEMCAL) {
103	fHistJetsPtClus = new TH1F("fHistJetsPtClus", "fHistJetsPtClus", fNbins, fMinPt, fMaxPt);
104	fHistJetsPtClus->GetXaxis()->SetTitle("p_{T} [GeV/c]");
105	fHistJetsPtClus->GetYaxis()->SetTitle("counts");
106	fOutput->Add(fHistJetsPtClus);
107	}
108
109	fHistJetsPtTrack = new TH1F("fHistJetsPtTrack", "fHistJetsPtTrack", fNbins, fMinPt, fMaxPt);
110	fHistJetsPtTrack->GetXaxis()->SetTitle("p_{T} [GeV/c]");
111	fHistJetsPtTrack->GetYaxis()->SetTitle("counts");
112	fOutput->Add(fHistJetsPtTrack);
113
114	fHistJetsPtNonBias = new TH1F("fHistJetsPtNonBias", "fHistJetsPtNonBias", fNbins, fMinPt, fMaxPt);
115	fHistJetsPtNonBias->GetXaxis()->SetTitle("p_{T} [GeV/c]");
116	fHistJetsPtNonBias->GetYaxis()->SetTitle("counts");
117	fOutput->Add(fHistJetsPtNonBias);
118
119	fHistLeadingJetPt = new TH1F("fHistLeadingJetPt", "fHistLeadingJetPt", fNbins, fMinPt, fMaxPt);
120	fHistLeadingJetPt->GetXaxis()->SetTitle("p_{T} [GeV]");
121	fHistLeadingJetPt->GetYaxis()->SetTitle("counts");
122	fOutput->Add(fHistLeadingJetPt);
123
124	fHistJetsZvsPt = new TH2F("fHistJetsZvsPt", "fHistJetsZvsPt", fNbins, 0, 1.2, fNbins, fMinPt, fMaxPt);
125	fHistJetsZvsPt->GetXaxis()->SetTitle("Z");
126	fHistJetsZvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
127	fOutput->Add(fHistJetsZvsPt);
128
129	if (fAnaType == kEMCAL) {
130	fHistJetsNEFvsPt = new TH2F("fHistJetsNEFvsPt", "fHistJetsNEFvsPt", fNbins, 0, 1.2, fNbins, fMinPt, fMaxPt);
131	fHistJetsNEFvsPt->GetXaxis()->SetTitle("NEF");
132	fHistJetsNEFvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
133	fOutput->Add(fHistJetsNEFvsPt);
134	}
135	}
136
137	//________________________________________________________________________
138	void AliJetResponseMaker::FillHistograms()
139	{
140	// Fill histograms.
141
142	Int_t maxJetIndex = -1;
143	Int_t maxMCJetIndex = -1;
144
145	DoJetLoop(maxJetIndex, fJets, fTracks, fCaloClusters);
146	if (maxJetIndex < 0)
147	return;
148	AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(fJets->At(maxJetIndex));
149	if (!jet)
150	return;
151
152	DoJetLoop(maxMCJetIndex, fMCJets, fMCParts);
153	if (maxMCJetIndex < 0)
154	return;
155	AliEmcalJet* mcjet = dynamic_cast<AliEmcalJet*>(fMCJets->At(maxMCJetIndex));
156	if (!mcjet)
157	return;
158
159	fHistLeadingJetPt->Fill(jet->Pt());
160	fHistMCLeadingJetPt->Fill(mcjet->Pt());
161	}
162
163	//________________________________________________________________________
164	void AliJetResponseMaker::DoJetLoop(Int_t &maxJetIndex, TClonesArray jets, TClonesArray tracks, TClonesArray *clusters)
165	{
166	// Do the jet loop.
167
168	if (!fJets)
169	return;
170
171	Int_t njets = jets->GetEntriesFast();
172
173	Float_t maxJetPt = 0;
174	for (Int_t ij = 0; ij < njets; ij++) {
175
176	AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(jets->At(ij));
177
178	if (!jet) {
179	AliError(Form("Could not receive jet %d", ij));
180	continue;
181	}
182
183	if (!AcceptJet(jet))
184	continue;
185
186	fHistJetsPtNonBias->Fill(jet->Pt());
187
188	if (jet->MaxTrackPt() > fPtBiasJetTrack)
189	fHistJetsPtTrack->Fill(jet->Pt());
190
191	if (fAnaType == kEMCAL && jet->MaxClusterPt() > fPtBiasJetClus)
192	fHistJetsPtClus->Fill(jet->Pt());
193
194	if (jet->MaxTrackPt() < fPtBiasJetTrack && (fAnaType == kTPC \|\| jet->MaxClusterPt() < fPtBiasJetClus))
195	continue;
196
197	if (maxJetPt < jet->Pt()) {
198	maxJetPt = jet->Pt();
199	maxJetIndex = ij;
200	}
201
202	fHistJetsPt->Fill(jet->Pt());
203
204	fHistJetPhiEta->Fill(jet->Eta(), jet->Phi());
205
206	if (fAnaType == kEMCAL)
207	fHistJetsNEFvsPt->Fill(jet->NEF(), jet->Pt());
208
209	for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
210	AliVParticle *track = jet->TrackAt(it, tracks);
211	if (track)
212	fHistJetsZvsPt->Fill(track->Pt() / jet->Pt(), jet->Pt());
213	}
214
215	if (fAnaType != kEMCAL \|\| !clusters)
216	continue;
217
218	for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
219	AliVCluster *cluster = jet->ClusterAt(ic, clusters);
220
221	if (cluster) {
222	TLorentzVector nPart;
223	cluster->GetMomentum(nPart, fVertex);
224	fHistJetsZvsPt->Fill(nPart.Et() / jet->Pt(), jet->Pt());
225	}
226	}
227	}
228	}
229
230	//________________________________________________________________________
231	void AliJetResponseMaker::RetrieveEventObjects()
232	{
233	// Retrieve event objects.
234
235	AliAnalysisTaskEmcal::RetrieveEventObjects();
236
237	if (strcmp(fMCJetsName,"")) {
238	fMCJets = dynamic_cast<TClonesArray*>(MCEvent()->FindListObject(fMCJetsName));
239	if (!fMCJets) {
240	AliWarning(Form("Could not retrieve MC jets %s!", fMCJetsName.Data()));
241	}
242	}
243	}
244
245	//________________________________________________________________________
246	void AliJetResponseMaker::Terminate(Option_t *)
247	{
248	// Called once at the end of the analysis.
249	}