3 // Emcal jet response matrix maker task.
7 #include "AliJetResponseMaker.h"
10 #include <TClonesArray.h>
14 #include <TLorentzVector.h>
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 "AliGenPythiaEventHeader.h"
23 #include "AliMCEvent.h"
25 ClassImp(AliJetResponseMaker)
27 //________________________________________________________________________
28 AliJetResponseMaker::AliJetResponseMaker() :
29 AliAnalysisTaskEmcalJet("AliJetResponseMaker", kTRUE),
30 fMCTracksName("MCParticles"),
31 fMCJetsName("MCJets"),
42 fHistAcceptedEvents(0),
46 fHistMCJetsNEFvsPt(0),
52 fHistClosestDistance(0),
53 fHistClosestDeltaPhi(0),
54 fHistClosestDeltaEta(0),
55 fHistClosestDeltaPt(0),
56 fHistNonMatchedMCJetPt(0),
57 fHistNonMatchedJetPt(0),
58 fHistPartvsDetecPt(0),
61 // Default constructor.
64 //________________________________________________________________________
65 AliJetResponseMaker::AliJetResponseMaker(const char *name) :
66 AliAnalysisTaskEmcalJet(name, kTRUE),
67 fMCTracksName("MCParticles"),
68 fMCJetsName("MCJets"),
79 fHistAcceptedEvents(0),
83 fHistMCJetsNEFvsPt(0),
89 fHistClosestDistance(0),
90 fHistClosestDeltaPhi(0),
91 fHistClosestDeltaEta(0),
92 fHistClosestDeltaPt(0),
93 fHistNonMatchedMCJetPt(0),
94 fHistNonMatchedJetPt(0),
95 fHistPartvsDetecPt(0),
98 // Standard constructor.
101 //________________________________________________________________________
102 AliJetResponseMaker::~AliJetResponseMaker()
107 //________________________________________________________________________
108 void AliJetResponseMaker::UserCreateOutputObjects()
110 // Create user objects.
113 fOutput = new TList();
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};
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);
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);
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);
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]));
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);
145 fHistJetsPt = new TH1F("fHistJetsPt", "fHistJetsPt", fNbins, fMinBinPt, fMaxBinPt);
146 fHistJetsPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
147 fHistJetsPt->GetYaxis()->SetTitle("counts");
148 fOutput->Add(fHistJetsPt);
150 fHistJetsZvsPt = new TH2F("fHistJetsZvsPt", "fHistJetsZvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
151 fHistJetsZvsPt->GetXaxis()->SetTitle("Z");
152 fHistJetsZvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
153 fOutput->Add(fHistJetsZvsPt);
155 if (fAnaType == kEMCAL) {
156 fHistJetsNEFvsPt = new TH2F("fHistJetsNEFvsPt", "fHistJetsNEFvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
157 fHistJetsNEFvsPt->GetXaxis()->SetTitle("NEF");
158 fHistJetsNEFvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
159 fOutput->Add(fHistJetsNEFvsPt);
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);
167 fHistMCJetsPt = new TH1F("fHistMCJetsPt", "fHistMCJetsPt", fNbins, fMinBinPt, fMaxBinPt);
168 fHistMCJetsPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
169 fHistMCJetsPt->GetYaxis()->SetTitle("counts");
170 fOutput->Add(fHistMCJetsPt);
172 fHistMCJetsZvsPt = new TH2F("fHistMCJetsZvsPt", "fHistMCJetsZvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
173 fHistMCJetsZvsPt->GetXaxis()->SetTitle("Z");
174 fHistMCJetsZvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
175 fOutput->Add(fHistMCJetsZvsPt);
177 if (fAnaType == kEMCAL) {
178 fHistMCJetsNEFvsPt = new TH2F("fHistMCJetsNEFvsPt", "fHistMCJetsNEFvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
179 fHistMCJetsNEFvsPt->GetXaxis()->SetTitle("NEF");
180 fHistMCJetsNEFvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
181 fOutput->Add(fHistMCJetsNEFvsPt);
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);
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);
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);
199 fHistClosestDeltaPt = new TH1F("fHistClosestDeltaPt", "fHistClosestDeltaPt", fNbins, -fMaxBinPt / 2, fMaxBinPt / 2);
200 fHistClosestDeltaPt->GetXaxis()->SetTitle("#Delta p_{T}");
201 fHistClosestDeltaPt->GetYaxis()->SetTitle("counts");
202 fOutput->Add(fHistClosestDeltaPt);
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);
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);
214 fHistPartvsDetecPt = new TH2F("fHistPartvsDetecPt", "fHistPartvsDetecPt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
215 fHistPartvsDetecPt->GetXaxis()->SetTitle("p_{T}^{rec}");
216 fHistPartvsDetecPt->GetYaxis()->SetTitle("p_{T}^{gen}");
217 fOutput->Add(fHistPartvsDetecPt);
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);
224 PostData(1, fOutput); // Post data for ALL output slots > 0 here, to get at least an empty histogram
227 //________________________________________________________________________
228 void AliJetResponseMaker::DoJetLoop(TClonesArray *jets1, TClonesArray *jets2, Bool_t mc)
232 Int_t nJets1 = jets1->GetEntriesFast();
233 Int_t nJets2 = jets2->GetEntriesFast();
235 for (Int_t i = 0; i < nJets1; i++) {
237 AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(jets1->At(i));
240 AliError(Form("Could not receive jet %d", i));
244 if (!AcceptJet(jet1, kTRUE, mc))
247 for (Int_t j = 0; j < nJets2; j++) {
249 AliEmcalJet* jet2 = static_cast<AliEmcalJet*>(jets2->At(j));
252 AliError(Form("Could not receive jet %d", j));
256 if (!AcceptJet(jet2, kTRUE, !mc))
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);
263 if (d < jet1->ClosestJetDistance()) {
264 jet1->SetSecondClosestJet(jet1->ClosestJet(), jet1->ClosestJetDistance());
265 jet1->SetClosestJet(jet2, d);
267 else if (d < jet1->SecondClosestJetDistance()) {
268 jet1->SetSecondClosestJet(jet2, d);
274 //________________________________________________________________________
275 void AliJetResponseMaker::ExecOnce()
279 if (!fMCJetsName.IsNull() && !fMCJets) {
280 fMCJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCJetsName));
282 AliError(Form("%s: Could not retrieve mc jets %s!", GetName(), fMCJetsName.Data()));
285 else if (!fMCJets->GetClass()->GetBaseClass("AliEmcalJet")) {
286 AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fMCJetsName.Data()));
292 if (!fMCTracksName.IsNull() && !fMCTracks) {
293 fMCTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCTracksName));
295 AliError(Form("%s: Could not retrieve mc tracks %s!", GetName(), fMCTracksName.Data()));
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()));
308 AliAnalysisTaskEmcalJet::ExecOnce();
311 //________________________________________________________________________
312 Bool_t AliJetResponseMaker::FillHistograms()
316 const Int_t nMCJets = fMCJets->GetEntriesFast();
318 for (Int_t i = 0; i < nMCJets; i++) {
320 AliEmcalJet* jet = static_cast<AliEmcalJet*>(fMCJets->At(i));
323 AliError(Form("Could not receive jet %d", i));
327 if (!AcceptJet(jet, kTRUE, kFALSE))
330 if (jet->Pt() > fMaxBinPt)
333 if (jet->ClosestJet() && jet->ClosestJet()->ClosestJet() == jet &&
334 jet->ClosestJetDistance() < fMaxDistance) { // Matched jet found
335 jet->SetMatchedToClosest();
336 jet->ClosestJet()->SetMatchedToClosest();
337 if (jet->MatchedJet()->Pt() > fMaxBinPt) {
338 fHistMissedMCJets->Fill(jet->Pt(), fEventWeight);
341 fHistClosestDistance->Fill(jet->ClosestJetDistance(), fEventWeight);
343 Double_t deta = jet->MatchedJet()->Eta() - jet->Eta();
344 fHistClosestDeltaEta->Fill(deta, fEventWeight);
346 Double_t dphi = jet->MatchedJet()->Phi() - jet->Phi();
347 fHistClosestDeltaPhi->Fill(dphi, fEventWeight);
349 Double_t dpt = jet->MatchedJet()->Pt() - jet->Pt();
350 fHistClosestDeltaPt->Fill(dpt, fEventWeight);
352 fHistPartvsDetecPt->Fill(jet->MatchedJet()->Pt(), jet->Pt(), fEventWeight);
356 fHistNonMatchedMCJetPt->Fill(jet->Pt(), fEventWeight);
357 fHistMissedMCJets->Fill(jet->Pt(), fEventWeight);
360 fHistMCJetsPt->Fill(jet->Pt(), fEventWeight);
362 fHistMCJetPhiEta->Fill(jet->Eta(), jet->Phi(), fEventWeight);
364 if (fAnaType == kEMCAL)
365 fHistMCJetsNEFvsPt->Fill(jet->NEF(), jet->Pt(), fEventWeight);
367 for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
368 AliVParticle *track = jet->TrackAt(it, fMCTracks);
370 fHistMCJetsZvsPt->Fill(track->Pt() / jet->Pt(), jet->Pt(), fEventWeight);
374 const Int_t nJets = fJets->GetEntriesFast();
376 for (Int_t i = 0; i < nJets; i++) {
378 AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(i));
381 AliError(Form("Could not receive mc jet %d", i));
388 if (!jet->MatchedJet()) {
389 fHistNonMatchedJetPt->Fill(jet->Pt(), fEventWeight);
392 fHistJetsPt->Fill(jet->Pt(), fEventWeight);
394 fHistJetPhiEta->Fill(jet->Eta(), jet->Phi(), fEventWeight);
396 if (fAnaType == kEMCAL)
397 fHistJetsNEFvsPt->Fill(jet->NEF(), jet->Pt(), fEventWeight);
399 for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
400 AliVParticle *track = jet->TrackAt(it, fTracks);
402 fHistJetsZvsPt->Fill(track->Pt() / jet->Pt(), jet->Pt(), fEventWeight);
405 for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
406 AliVCluster *cluster = jet->ClusterAt(ic, fCaloClusters);
409 cluster->GetMomentum(nP, fVertex);
410 fHistJetsZvsPt->Fill(nP.Pt() / jet->Pt(), jet->Pt(), fEventWeight);
418 //________________________________________________________________________
419 Bool_t AliJetResponseMaker::RetrieveEventObjects()
421 // Retrieve event objects.
423 if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
426 fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
432 fEventWeight = fPythiaHeader->EventWeight();
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};
439 Double_t pthard = fPythiaHeader->GetPtHard();
441 for (fPtHardBin = 0; fPtHardBin < 11; fPtHardBin++) {
442 if (ptHardLo[fPtHardBin] < pthard && ptHardHi[fPtHardBin] > pthard)
446 fNTrials = fPythiaHeader->Trials();
451 //________________________________________________________________________
452 Bool_t AliJetResponseMaker::Run()
454 // Find the closest jets
456 fHistEvents->SetBinContent(fPtHardBin, fHistEvents->GetBinContent(fPtHardBin) + 1);
458 if (TMath::Abs(fVertex[2]) > fVertexCut)
461 fHistAcceptedEvents->SetBinContent(fPtHardBin, fHistEvents->GetBinContent(fPtHardBin) + 1);
462 fHistNTrials->SetBinContent(fPtHardBin, fHistNTrials->GetBinContent(fPtHardBin) + fNTrials);
464 DoJetLoop(fJets, fMCJets, kFALSE);
465 DoJetLoop(fMCJets, fJets, kTRUE);