9 #include <TClonesArray.h>
13 #include <TLorentzVector.h>
15 #include <TParameter.h>
17 #include "AliAnalysisManager.h"
18 #include "AliCentrality.h"
19 #include "AliVCluster.h"
20 #include "AliVParticle.h"
21 #include "AliVTrack.h"
22 #include "AliEmcalJet.h"
23 #include "AliVEventHandler.h"
26 #include "AliAnalysisTaskSAJF.h"
28 ClassImp(AliAnalysisTaskSAJF)
30 //________________________________________________________________________
31 AliAnalysisTaskSAJF::AliAnalysisTaskSAJF() :
32 AliAnalysisTaskEmcalJet("AliAnalysisTaskSAJF"),
34 fEmbJetsName("EmbJets"),
35 fRandTracksName("TracksRandomized"),
36 fRandCaloName("CaloClustersRandomized"),
38 fSkipEventsWithNoJets(kTRUE),
44 fHistRejectedEvents(0),
45 fHistRhoVSleadJetPt(0),
47 fHistRCPtExLJVSDPhiLJ(0),
50 fHistEmbPartPhiEta(0),
54 // Default constructor.
56 for (Int_t i = 0; i < 4; i++) {
57 fHistJetPhiEta[i] = 0;
59 fHistJetsPtArea[i] = 0;
60 fHistJetsPtNonBias[i] = 0;
61 fHistJetsPtAreaNonBias[i] = 0;
62 fHistLeadingJetPt[i] = 0;
63 fHist2LeadingJetPt[i] = 0;
64 fHistJetsNEFvsPt[i] = 0;
65 fHistJetsZvsPt[i] = 0;
67 fHistCorrJetsPt[i] = 0;
68 fHistCorrJetsPtArea[i] = 0;
69 fHistCorrJetsPtNonBias[i] = 0;
70 fHistCorrJetsPtAreaNonBias[i] = 0;
71 fHistCorrLeadingJetPt[i] = 0;
75 fHistDeltaPtRC[i] = 0;
76 fHistDeltaPtRCExLJ[i] = 0;
77 fHistDeltaPtRCRand[i] = 0;
80 fHistDeltaPtEmb[i] = 0;
84 //________________________________________________________________________
85 AliAnalysisTaskSAJF::AliAnalysisTaskSAJF(const char *name) :
86 AliAnalysisTaskEmcalJet(name),
88 fEmbJetsName("EmbJets"),
89 fRandTracksName("TracksRandomized"),
90 fRandCaloName("CaloClustersRandomized"),
92 fSkipEventsWithNoJets(kTRUE),
98 fHistRejectedEvents(0),
99 fHistRhoVSleadJetPt(0),
101 fHistRCPtExLJVSDPhiLJ(0),
103 fHistEmbJetPhiEta(0),
104 fHistEmbPartPhiEta(0),
107 // Standard constructor.
109 for (Int_t i = 0; i < 4; i++) {
110 fHistJetPhiEta[i] = 0;
112 fHistJetsPtArea[i] = 0;
113 fHistJetsPtNonBias[i] = 0;
114 fHistJetsPtAreaNonBias[i] = 0;
115 fHistLeadingJetPt[i] = 0;
116 fHist2LeadingJetPt[i] = 0;
117 fHistJetsNEFvsPt[i] = 0;
118 fHistJetsZvsPt[i] = 0;
120 fHistCorrJetsPt[i] = 0;
121 fHistCorrJetsPtArea[i] = 0;
122 fHistCorrJetsPtNonBias[i] = 0;
123 fHistCorrJetsPtAreaNonBias[i] = 0;
124 fHistCorrLeadingJetPt[i] = 0;
126 fHistRCPtExLJ[i] = 0;
127 fHistRCPtRand[i] = 0;
128 fHistDeltaPtRC[i] = 0;
129 fHistDeltaPtRCExLJ[i] = 0;
130 fHistDeltaPtRCRand[i] = 0;
133 fHistDeltaPtEmb[i] = 0;
137 //________________________________________________________________________
138 AliAnalysisTaskSAJF::~AliAnalysisTaskSAJF()
143 //________________________________________________________________________
144 void AliAnalysisTaskSAJF::UserCreateOutputObjects()
146 // Create user output.
148 AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
150 const Float_t binWidth = (fMaxBinPt - fMinBinPt) / fNbins;
153 fOutput = new TList();
156 fHistCentrality = new TH1F("fHistCentrality","fHistCentrality", fNbins, 0, 100);
157 fHistCentrality->GetXaxis()->SetTitle("Centrality (%)");
158 fHistCentrality->GetYaxis()->SetTitle("counts");
159 fOutput->Add(fHistCentrality);
161 fHistRejectedEvents = new TH1F("fHistRejectedEvents","fHistRejectedEvents", 6, 0, 6);
162 fHistRejectedEvents->GetXaxis()->SetTitle("Reason");
163 fHistRejectedEvents->GetYaxis()->SetTitle("counts");
164 fHistRejectedEvents->GetXaxis()->SetBinLabel(1, "Rho <= 0");
165 fHistRejectedEvents->GetXaxis()->SetBinLabel(2, "Max Jet <= 0");
166 fHistRejectedEvents->GetXaxis()->SetBinLabel(3, "Max Jet not found");
167 fHistRejectedEvents->GetXaxis()->SetBinLabel(4, "No jets");
168 fOutput->Add(fHistRejectedEvents);
170 fHistRhoVSleadJetPt = new TH2F("fHistRhoVSleadJetPt","fHistRhoVSleadJetPt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
171 fHistRhoVSleadJetPt->GetXaxis()->SetTitle("#rho * area [GeV/c]");
172 fHistRhoVSleadJetPt->GetYaxis()->SetTitle("Leading jet p_{T} [GeV/c]");
173 fOutput->Add(fHistRhoVSleadJetPt);
175 fHistRCPhiEta = new TH2F("fHistRCPhiEta","Phi-Eta distribution of rigid cones", 40, -2, 2, 64, 0, 6.4);
176 fHistRCPhiEta->GetXaxis()->SetTitle("#eta");
177 fHistRCPhiEta->GetYaxis()->SetTitle("#phi");
178 fOutput->Add(fHistRCPhiEta);
180 fHistRCPtExLJVSDPhiLJ = new TH2F("fHistRCPtExLJVSDPhiLJ","fHistRCPtExLJVSDPhiLJ", fNbins, fMinBinPt, fMaxBinPt, 128, -1.6, 4.8);
181 fHistRCPtExLJVSDPhiLJ->GetXaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
182 fHistRCPtExLJVSDPhiLJ->GetYaxis()->SetTitle("#Delta#phi");
183 fOutput->Add(fHistRCPtExLJVSDPhiLJ);
185 fHistRhoVSRCPt = new TH2F("fHistRhoVSRCPt","fHistRhoVSRCPt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
186 fHistRhoVSRCPt->GetXaxis()->SetTitle("#rho * area [GeV/c]");
187 fHistRhoVSRCPt->GetYaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
188 fOutput->Add(fHistRhoVSRCPt);
190 fHistEmbJetPhiEta = new TH2F("fHistEmbJetPhiEta","Phi-Eta distribution of embedded jets", 40, -2, 2, 64, 0, 6.4);
191 fHistEmbJetPhiEta->GetXaxis()->SetTitle("#eta");
192 fHistEmbJetPhiEta->GetYaxis()->SetTitle("#phi");
193 fOutput->Add(fHistEmbJetPhiEta);
195 fHistEmbPartPhiEta = new TH2F("fHistEmbPartPhiEta","Phi-Eta distribution of embedded particles", 40, -2, 2, 64, 0, 6.4);
196 fHistEmbPartPhiEta->GetXaxis()->SetTitle("#eta");
197 fHistEmbPartPhiEta->GetYaxis()->SetTitle("#phi");
198 fOutput->Add(fHistEmbPartPhiEta);
200 fHistRhoVSEmbBkg = new TH2F("fHistRhoVSEmbBkg","fHistRhoVSEmbBkg", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
201 fHistRhoVSEmbBkg->GetXaxis()->SetTitle("rho * area [GeV/c]");
202 fHistRhoVSEmbBkg->GetYaxis()->SetTitle("background of embedded track [GeV/c]");
203 fOutput->Add(fHistRhoVSEmbBkg);
207 for (Int_t i = 0; i < 4; i++) {
208 histname = "fHistJetPhiEta_";
210 fHistJetPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 40, -2, 2, 64, 0, 6.4);
211 fHistJetPhiEta[i]->GetXaxis()->SetTitle("#eta");
212 fHistJetPhiEta[i]->GetYaxis()->SetTitle("#phi");
213 fOutput->Add(fHistJetPhiEta[i]);
215 histname = "fHistJetsPt_";
217 fHistJetsPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
218 fHistJetsPt[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
219 fHistJetsPt[i]->GetYaxis()->SetTitle("counts");
220 fOutput->Add(fHistJetsPt[i]);
222 histname = "fHistJetsPtArea_";
224 fHistJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, 20, 0, fJetRadius * fJetRadius * TMath::Pi() * 1.5);
225 fHistJetsPtArea[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
226 fHistJetsPtArea[i]->GetYaxis()->SetTitle("area");
227 fOutput->Add(fHistJetsPtArea[i]);
229 histname = "fHistJetsPtNonBias_";
231 fHistJetsPtNonBias[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
232 fHistJetsPtNonBias[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
233 fHistJetsPtNonBias[i]->GetYaxis()->SetTitle("counts");
234 fOutput->Add(fHistJetsPtNonBias[i]);
236 histname = "fHistJetsPtAreaNonBias_";
238 fHistJetsPtAreaNonBias[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, 20, 0, fJetRadius * fJetRadius * TMath::Pi() * 1.5);
239 fHistJetsPtAreaNonBias[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
240 fHistJetsPtAreaNonBias[i]->GetYaxis()->SetTitle("area");
241 fOutput->Add(fHistJetsPtAreaNonBias[i]);
243 histname = "fHistLeadingJetPt_";
245 fHistLeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
246 fHistLeadingJetPt[i]->GetXaxis()->SetTitle("p_{T} [GeV]");
247 fHistLeadingJetPt[i]->GetYaxis()->SetTitle("counts");
248 fOutput->Add(fHistLeadingJetPt[i]);
250 histname = "fHist2LeadingJetPt_";
252 fHist2LeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
253 fHist2LeadingJetPt[i]->GetXaxis()->SetTitle("p_{T} [GeV]");
254 fHist2LeadingJetPt[i]->GetYaxis()->SetTitle("counts");
255 fOutput->Add(fHist2LeadingJetPt[i]);
257 histname = "fHistJetsZvsPt_";
259 fHistJetsZvsPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
260 fHistJetsZvsPt[i]->GetXaxis()->SetTitle("Z");
261 fHistJetsZvsPt[i]->GetYaxis()->SetTitle("p_{T} [GeV/c]");
262 fOutput->Add(fHistJetsZvsPt[i]);
264 if (fAnaType == kEMCAL) {
265 histname = "fHistJetsNEFvsPt_";
267 fHistJetsNEFvsPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
268 fHistJetsNEFvsPt[i]->GetXaxis()->SetTitle("NEF");
269 fHistJetsNEFvsPt[i]->GetYaxis()->SetTitle("p_{T} [GeV/c]");
270 fOutput->Add(fHistJetsNEFvsPt[i]);
273 histname = "fHistRho_";
275 fHistRho[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
276 fHistRho[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
277 fHistRho[i]->GetYaxis()->SetTitle("counts");
278 fOutput->Add(fHistRho[i]);
280 histname = "fHistCorrJetsPt_";
282 fHistCorrJetsPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
283 fHistCorrJetsPt[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
284 fHistCorrJetsPt[i]->GetYaxis()->SetTitle("counts");
285 fOutput->Add(fHistCorrJetsPt[i]);
287 histname = "fHistCorrJetsPtArea_";
289 fHistCorrJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, 20, 0, fJetRadius * fJetRadius * TMath::Pi() * 1.5);
290 fHistCorrJetsPtArea[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
291 fHistCorrJetsPtArea[i]->GetYaxis()->SetTitle("area");
292 fOutput->Add(fHistCorrJetsPtArea[i]);
294 histname = "fHistCorrJetsPtNonBias_";
296 fHistCorrJetsPtNonBias[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
297 fHistCorrJetsPtNonBias[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
298 fHistCorrJetsPtNonBias[i]->GetYaxis()->SetTitle("counts");
299 fOutput->Add(fHistCorrJetsPtNonBias[i]);
301 histname = "fHistCorrJetsPtAreaNonBias_";
303 fHistCorrJetsPtAreaNonBias[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, 20, 0, fJetRadius * fJetRadius * TMath::Pi() * 1.5);
304 fHistCorrJetsPtAreaNonBias[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
305 fHistCorrJetsPtAreaNonBias[i]->GetYaxis()->SetTitle("area");
306 fOutput->Add(fHistCorrJetsPtAreaNonBias[i]);
308 histname = "fHistCorrLeadingJetPt_";
310 fHistCorrLeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
311 fHistCorrLeadingJetPt[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
312 fHistCorrLeadingJetPt[i]->GetYaxis()->SetTitle("counts");
313 fOutput->Add(fHistCorrLeadingJetPt[i]);
315 histname = "fHistRCPt_";
317 fHistRCPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
318 fHistRCPt[i]->GetXaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
319 fHistRCPt[i]->GetYaxis()->SetTitle("counts");
320 fOutput->Add(fHistRCPt[i]);
322 histname = "fHistRCPtExLJ_";
324 fHistRCPtExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
325 fHistRCPtExLJ[i]->GetXaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
326 fHistRCPtExLJ[i]->GetYaxis()->SetTitle("counts");
327 fOutput->Add(fHistRCPtExLJ[i]);
329 histname = "fHistRCPtRand_";
331 fHistRCPtRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
332 fHistRCPtRand[i]->GetXaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
333 fHistRCPtRand[i]->GetYaxis()->SetTitle("counts");
334 fOutput->Add(fHistRCPtRand[i]);
336 histname = "fHistDeltaPtRC_";
338 fHistDeltaPtRC[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt - fMaxBinPt / 2 + binWidth / 2, fMinBinPt + fMaxBinPt / 2 + binWidth / 2);
339 fHistDeltaPtRC[i]->GetXaxis()->SetTitle("#deltap_{T} [GeV/c]");
340 fHistDeltaPtRC[i]->GetYaxis()->SetTitle("counts");
341 fOutput->Add(fHistDeltaPtRC[i]);
343 histname = "fHistDeltaPtRCExLJ_";
345 fHistDeltaPtRCExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt - fMaxBinPt / 2 + binWidth / 2, fMinBinPt + fMaxBinPt / 2 + binWidth / 2);
346 fHistDeltaPtRCExLJ[i]->GetXaxis()->SetTitle("#deltap_{T} [GeV/c]");
347 fHistDeltaPtRCExLJ[i]->GetYaxis()->SetTitle("counts");
348 fOutput->Add(fHistDeltaPtRCExLJ[i]);
350 histname = "fHistDeltaPtRCRand_";
352 fHistDeltaPtRCRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt - fMaxBinPt / 2 + binWidth / 2, fMinBinPt + fMaxBinPt / 2 + binWidth / 2);
353 fHistDeltaPtRCRand[i]->GetXaxis()->SetTitle("#deltap_{T} [GeV/c]");
354 fHistDeltaPtRCRand[i]->GetYaxis()->SetTitle("counts");
355 fOutput->Add(fHistDeltaPtRCRand[i]);
357 histname = "fHistEmbJets_";
359 fHistEmbJets[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
360 fHistEmbJets[i]->GetXaxis()->SetTitle("embedded jet p_{T} [GeV/c]");
361 fHistEmbJets[i]->GetYaxis()->SetTitle("counts");
362 fOutput->Add(fHistEmbJets[i]);
364 histname = "fHistEmbPart_";
366 fHistEmbPart[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
367 fHistEmbPart[i]->GetXaxis()->SetTitle("embedded particle p_{T} [GeV/c]");
368 fHistEmbPart[i]->GetYaxis()->SetTitle("counts");
369 fOutput->Add(fHistEmbPart[i]);
371 histname = "fHistDeltaPtEmb_";
373 fHistDeltaPtEmb[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt - fMaxBinPt / 2 + binWidth / 2, fMinBinPt + fMaxBinPt / 2 + binWidth / 2);
374 fHistDeltaPtEmb[i]->GetXaxis()->SetTitle("#deltap_{T} [GeV/c]");
375 fHistDeltaPtEmb[i]->GetYaxis()->SetTitle("counts");
376 fOutput->Add(fHistDeltaPtEmb[i]);
379 PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
382 //________________________________________________________________________
383 Bool_t AliAnalysisTaskSAJF::RetrieveEventObjects()
385 // Retrieve event objects.
387 if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
392 if (strcmp(fRhoName,"")) {
393 TParameter<Double_t> *rho = dynamic_cast<TParameter<Double_t>*>(InputEvent()->FindListObject(fRhoName));
396 fRho = rho->GetVal();
399 AliWarning(Form("Could not retrieve rho %s!", fRhoName.Data()));
403 if (strcmp(fEmbJetsName,"")) {
404 fEmbJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbJetsName));
406 AliWarning(Form("Could not retrieve emb jets %s!", fEmbJetsName.Data()));
410 if (strcmp(fRandCaloName,"") && fAnaType == kEMCAL) {
411 fRandCaloClusters = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fRandCaloName));
412 if (!fRandCaloClusters) {
413 AliWarning(Form("Could not retrieve randomized clusters %s!", fRandCaloName.Data()));
417 if (strcmp(fRandTracksName,"")) {
418 fRandTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fRandTracksName));
420 AliWarning(Form("Could not retrieve randomized tracks %s!", fRandTracksName.Data()));
427 //________________________________________________________________________
428 Bool_t AliAnalysisTaskSAJF::FillHistograms()
433 fHistRejectedEvents->Fill("Rho <= 0", 1);
437 Int_t maxJetIndex = -1;
438 Int_t max2JetIndex = -1;
441 // General histograms
442 // _________________________________
444 GetLeadingJets(maxJetIndex, max2JetIndex);
446 if (fSkipEventsWithNoJets && maxJetIndex < 0) {
447 fHistRejectedEvents->Fill("No jets", 1);
451 AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(fJets->At(maxJetIndex));
452 if (fSkipEventsWithNoJets && !jet) {
453 fHistRejectedEvents->Fill("Max Jet not found", 1);
457 Float_t maxJetCorrPt = 0;
460 maxJetCorrPt = jet->Pt() - fRho * jet->Area();
462 if (fSkipEventsWithNoJets && maxJetCorrPt <= 0) {
463 fHistRejectedEvents->Fill("Max Jet <= 0", 1);
467 fHistCentrality->Fill(fCent);
468 fHistRho[fCentBin]->Fill(fRho);
471 fHistLeadingJetPt[fCentBin]->Fill(jet->Pt());
472 fHistRhoVSleadJetPt->Fill(fRho * jet->Area(), jet->Pt());
473 fHistCorrLeadingJetPt[fCentBin]->Fill(maxJetCorrPt);
476 AliEmcalJet* jet2 = 0;
477 if (max2JetIndex >= 0)
478 jet2 = dynamic_cast<AliEmcalJet*>(fJets->At(max2JetIndex));
481 fHist2LeadingJetPt[fCentBin]->Fill(jet2->Pt());
487 // _________________________________
489 const Float_t rcArea = fJetRadius * fJetRadius * TMath::Pi();
491 // Simple random cones
495 GetRigidCone(RCpt, RCeta, RCphi, kFALSE, 0);
497 fHistRCPt[fCentBin]->Fill(RCpt / rcArea);
498 fHistDeltaPtRC[fCentBin]->Fill(RCpt - rcArea * fRho);
501 // Random cones far from leading jet
502 Float_t RCptExLJ = 0;
503 Float_t RCetaExLJ = 0;
504 Float_t RCphiExLJ = 0;
505 GetRigidCone(RCptExLJ, RCetaExLJ, RCphiExLJ, kFALSE, jet);
507 fHistRCPhiEta->Fill(RCetaExLJ, RCphiExLJ);
508 fHistRhoVSRCPt->Fill(fRho, RCptExLJ / rcArea);
510 Float_t dphi = RCphiExLJ - jet->Phi();
511 if (dphi > 4.8) dphi -= TMath::Pi() * 2;
512 if (dphi < -1.6) dphi += TMath::Pi() * 2;
513 fHistRCPtExLJVSDPhiLJ->Fill(RCptExLJ, dphi);
515 fHistRCPtExLJ[fCentBin]->Fill(RCptExLJ / rcArea);
516 fHistDeltaPtRCExLJ[fCentBin]->Fill(RCptExLJ - rcArea * fRho);
519 // Random cones with randomized particles
520 Float_t RCptRand = 0;
521 Float_t RCetaRand = 0;
522 Float_t RCphiRand = 0;
523 GetRigidCone(RCptRand, RCetaRand, RCphiRand, kTRUE, 0, fRandTracks, fRandCaloClusters);
525 fHistRCPtRand[fCentBin]->Fill(RCptRand / rcArea);
526 fHistDeltaPtRCRand[fCentBin]->Fill(RCptRand - rcArea * fRho);
531 // _________________________________
536 AliEmcalJet *embJet = 0;
537 TObject *maxPart = 0;
539 DoEmbJetLoop(embJet, maxPart);
543 fHistEmbJets[fCentBin]->Fill(embJet->Pt());
544 fHistEmbPart[fCentBin]->Fill(embJet->MCPt());
545 fHistEmbJetPhiEta->Fill(embJet->Eta(), embJet->Phi());
546 fHistDeltaPtEmb[fCentBin]->Fill(embJet->Pt() - embJet->Area() * fRho - embJet->MCPt());
547 fHistRhoVSEmbBkg->Fill(embJet->Area() * fRho, embJet->Pt() - embJet->MCPt());
549 AliVCluster *cluster = dynamic_cast<AliVCluster*>(maxPart);
552 cluster->GetPosition(pos);
553 TVector3 clusVec(pos);
555 fHistEmbPartPhiEta->Fill(clusVec.Eta(), clusVec.Phi());
558 AliVTrack *track = dynamic_cast<AliVTrack*>(maxPart);
560 fHistEmbPartPhiEta->Fill(track->Eta(), track->Phi());
563 AliWarning(Form("%s - Embedded particle type not found or not recognized (neither AliVCluster nor AliVTrack)!", GetName()));
570 AliWarning(Form("%s - Embedded jet not found in the event!", GetName()));
576 //________________________________________________________________________
577 void AliAnalysisTaskSAJF::GetLeadingJets(Int_t &maxJetIndex, Int_t &max2JetIndex)
579 // Get the leading jets.
584 const Int_t njets = fJets->GetEntriesFast();
586 Float_t maxJetPt = 0;
587 Float_t max2JetPt = 0;
588 for (Int_t ij = 0; ij < njets; ij++) {
590 AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(fJets->At(ij));
593 AliError(Form("Could not receive jet %d", ij));
600 Float_t corrPt = jet->Pt() - fRho * jet->Area();
602 if (maxJetIndex == -1 || maxJetPt < corrPt) {
603 max2JetPt = maxJetPt;
604 max2JetIndex = maxJetIndex;
608 else if (max2JetIndex == -1 || max2JetPt < corrPt) {
615 //________________________________________________________________________
616 void AliAnalysisTaskSAJF::DoJetLoop()
623 const Int_t njets = fJets->GetEntriesFast();
625 for (Int_t ij = 0; ij < njets; ij++) {
627 AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(fJets->At(ij));
630 AliError(Form("Could not receive jet %d", ij));
634 if (!AcceptJet(jet, kFALSE))
637 Float_t corrPt = jet->Pt() - fRho * jet->Area();
639 fHistJetsPtNonBias[fCentBin]->Fill(jet->Pt());
640 fHistJetsPtAreaNonBias[fCentBin]->Fill(jet->Pt(), jet->Area());
641 fHistCorrJetsPtNonBias[fCentBin]->Fill(corrPt);
642 fHistCorrJetsPtAreaNonBias[fCentBin]->Fill(corrPt, jet->Area());
644 if (!AcceptBiasJet(jet))
647 fHistJetsPt[fCentBin]->Fill(jet->Pt());
648 fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area());
649 fHistCorrJetsPt[fCentBin]->Fill(corrPt);
650 fHistCorrJetsPtArea[fCentBin]->Fill(corrPt, jet->Area());
652 fHistJetPhiEta[fCentBin]->Fill(jet->Eta(), jet->Phi());
654 if (fAnaType == kEMCAL)
655 fHistJetsNEFvsPt[fCentBin]->Fill(jet->NEF(), jet->Pt());
658 for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
659 AliVParticle *track = jet->TrackAt(it, fTracks);
661 fHistJetsZvsPt[fCentBin]->Fill(track->Pt() / jet->Pt(), jet->Pt());
666 for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
667 AliVCluster *cluster = jet->ClusterAt(ic, fCaloClusters);
670 TLorentzVector nPart;
671 cluster->GetMomentum(nPart, fVertex);
672 fHistJetsZvsPt[fCentBin]->Fill(nPart.Et() / jet->Pt(), jet->Pt());
679 //________________________________________________________________________
680 void AliAnalysisTaskSAJF::DoEmbJetLoop(AliEmcalJet* &embJet, TObject* &maxPart)
682 // Do the embedded jet loop.
687 TLorentzVector maxClusVect;
689 Int_t nembjets = fEmbJets->GetEntriesFast();
691 for (Int_t ij = 0; ij < nembjets; ij++) {
693 AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(fEmbJets->At(ij));
696 AliError(Form("Could not receive jet %d", ij));
706 AliVParticle *maxTrack = 0;
708 for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
709 AliVParticle *track = jet->TrackAt(it, fTracks);
714 if (track->GetLabel() != 100)
717 if (!maxTrack || track->Pt() > maxTrack->Pt())
721 AliVCluster *maxClus = 0;
723 for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
724 AliVCluster *cluster = jet->ClusterAt(ic, fCaloClusters);
729 if (cluster->Chi2() != 100)
732 TLorentzVector nPart;
733 cluster->GetMomentum(nPart, fVertex);
735 if (!maxClus || nPart.Et() > maxClusVect.Et()) {
736 new (&maxClusVect) TLorentzVector(nPart);
741 if ((maxClus && maxTrack && maxClusVect.Et() > maxTrack->Pt()) || (maxClus && !maxTrack)) {
750 return; // MC jets found, exit
754 //________________________________________________________________________
755 void AliAnalysisTaskSAJF::GetRigidCone(Float_t &pt, Float_t &eta, Float_t &phi, Bool_t acceptMC,
756 AliEmcalJet *jet, TClonesArray* tracks, TClonesArray* clusters) const
764 clusters = fCaloClusters;
766 if (!tracks && !clusters)
781 Float_t maxEta = fMaxEta;
782 Float_t minEta = fMinEta;
783 Float_t maxPhi = fMaxPhi;
784 Float_t minPhi = fMinPhi;
786 if (maxPhi > TMath::Pi() * 2) maxPhi = TMath::Pi() * 2;
787 if (minPhi < 0) minPhi = 0;
792 eta = gRandom->Rndm() * (maxEta - minEta) + minEta;
793 phi = gRandom->Rndm() * (maxPhi - minPhi) + minPhi;
794 dLJ = TMath::Sqrt((LJeta - eta) * (LJeta - eta) + (LJphi - phi) * (LJphi - phi));
796 } while (dLJ < fMinRC2LJ && repeats < 999);
801 if (fAnaType == kEMCAL && clusters) {
802 Int_t nclusters = clusters->GetEntriesFast();
803 for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
804 AliVCluster* cluster = dynamic_cast<AliVCluster*>(clusters->At(iClusters));
806 AliError(Form("Could not receive cluster %d", iClusters));
810 if (!(cluster->IsEMCAL())) continue;
812 if (!AcceptCluster(cluster, acceptMC)) continue;
814 TLorentzVector nPart;
815 cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
818 cluster->GetPosition(pos);
819 TVector3 clusVec(pos);
821 Float_t d = TMath::Sqrt((clusVec.Eta() - eta) * (clusVec.Eta() - eta) + (clusVec.Phi() - phi) * (clusVec.Phi() - phi));
829 Int_t ntracks = tracks->GetEntriesFast();
830 for(Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
831 AliVTrack* track = dynamic_cast<AliVTrack*>(tracks->At(iTracks));
833 AliError(Form("Could not retrieve track %d",iTracks));
837 if (!AcceptTrack(track, acceptMC)) continue;
839 Float_t tracketa = track->Eta();
840 Float_t trackphi = track->Phi();
842 if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi + 2 * TMath::Pi()))
843 trackphi += 2 * TMath::Pi();
844 if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi - 2 * TMath::Pi()))
845 trackphi -= 2 * TMath::Pi();
847 Float_t d = TMath::Sqrt((tracketa - eta) * (tracketa - eta) + (trackphi - phi) * (trackphi - phi));
854 //________________________________________________________________________
855 void AliAnalysisTaskSAJF::Init()
857 // Initialize the analysis.
859 AliAnalysisTaskEmcalJet::Init();
861 const Float_t semiDiag = TMath::Sqrt((fMaxPhi - fMinPhi) * (fMaxPhi - fMinPhi) +
862 (fMaxEta - fMinEta) * (fMaxEta - fMinEta)) / 2;
863 if (fMinRC2LJ > semiDiag * 0.5) {
864 AliWarning(Form("The parameter fMinRC2LJ = %f is too large for the considered acceptance. "
865 "Will use fMinRC2LJ = %f", fMinRC2LJ, semiDiag * 0.5));
866 fMinRC2LJ = semiDiag * 0.5;
870 //________________________________________________________________________
871 void AliAnalysisTaskSAJF::Terminate(Option_t *)
873 // Called once at the end of the analysis.