Changes by Salvatore
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliAnalysisTaskDeltaPt.cxx
1 // $Id$
2 //
3 // Jet deltaPt task.
4 //
5 // Author: S.Aiola
6
7 #include <TClonesArray.h>
8 #include <TH1F.h>
9 #include <TH2F.h>
10 #include <TList.h>
11 #include <TLorentzVector.h>
12 #include <TRandom3.h>
13
14 #include "AliVCluster.h"
15 #include "AliVParticle.h"
16 #include "AliVTrack.h"
17 #include "AliEmcalJet.h"
18 #include "AliRhoParameter.h"
19 #include "AliLog.h"
20
21 #include "AliAnalysisTaskDeltaPt.h"
22
23 ClassImp(AliAnalysisTaskDeltaPt)
24
25 //________________________________________________________________________
26 AliAnalysisTaskDeltaPt::AliAnalysisTaskDeltaPt() : 
27   AliAnalysisTaskEmcalJet("AliAnalysisTaskDeltaPt", kTRUE),
28   fMCAna(kFALSE),
29   fMinRC2LJ(-1),
30   fEmbJetsName(""),
31   fEmbTracksName(""),
32   fEmbCaloName(""),
33   fRandTracksName("TracksRandomized"),
34   fRandCaloName("CaloClustersRandomized"),
35   fRCperEvent(-1),
36   fEmbJets(0),
37   fEmbTracks(0),
38   fEmbCaloClusters(0),
39   fRandTracks(0),
40   fRandCaloClusters(0),
41   fEmbeddedClusterId(-1),
42   fEmbeddedTrackId(-1),
43   fHistRCPhiEta(0),
44   fHistRCPtExLJVSDPhiLJ(0),
45   fHistEmbJetPhiEta(0),
46   fHistEmbPartPhiEta(0)
47 {
48   // Default constructor.
49
50   for (Int_t i = 0; i < 4; i++) {
51     fHistRCPt[i] = 0;
52     fHistRCPtExLJ[i] = 0;
53     fHistRCPtRand[i] = 0; 
54     fHistRhoVSRCPt[i] = 0;
55     fHistDeltaPtRC[i] = 0;
56     fHistDeltaPtRCExLJ[i] = 0;
57     fHistDeltaPtRCRand[i] = 0;
58     fHistEmbNotFoundPhiEta[i] = 0;
59     fHistEmbJetsPtArea[i] = 0;
60     fHistEmbJetsCorrPtArea[i] = 0;
61     fHistEmbPartPt[i] = 0;
62     fHistDistEmbPartJetAxis[i] = 0;
63     fHistEmbBkgArea[i] = 0;
64     fHistRhoVSEmbBkg[i] = 0;
65     fHistDeltaPtEmbArea[i] = 0;
66   }
67
68   SetMakeGeneralHistograms(kTRUE);
69 }
70
71 //________________________________________________________________________
72 AliAnalysisTaskDeltaPt::AliAnalysisTaskDeltaPt(const char *name) : 
73   AliAnalysisTaskEmcalJet(name, kTRUE),
74   fMCAna(kFALSE),
75   fMinRC2LJ(-1),
76   fEmbJetsName(""),
77   fEmbTracksName(""),
78   fEmbCaloName(""),
79   fRandTracksName("TracksRandomized"),
80   fRandCaloName("CaloClustersRandomized"),
81   fRCperEvent(-1),
82   fEmbJets(0),
83   fEmbTracks(0),
84   fEmbCaloClusters(0),
85   fRandTracks(0),
86   fRandCaloClusters(0),
87   fEmbeddedClusterId(-1),
88   fEmbeddedTrackId(-1),
89   fHistRCPhiEta(0),
90   fHistRCPtExLJVSDPhiLJ(0),
91   fHistEmbJetPhiEta(0),
92   fHistEmbPartPhiEta(0)
93 {
94   // Standard constructor.
95
96   for (Int_t i = 0; i < 4; i++) {
97     fHistRCPt[i] = 0;
98     fHistRCPtExLJ[i] = 0;
99     fHistRCPtRand[i] = 0; 
100     fHistRhoVSRCPt[i] = 0;
101     fHistDeltaPtRC[i] = 0;
102     fHistDeltaPtRCExLJ[i] = 0;
103     fHistDeltaPtRCRand[i] = 0;
104     fHistEmbNotFoundPhiEta[i] = 0;
105     fHistEmbJetsPtArea[i] = 0;
106     fHistEmbJetsCorrPtArea[i] = 0;
107     fHistEmbPartPt[i] = 0;
108     fHistDistEmbPartJetAxis[i] = 0;
109     fHistEmbBkgArea[i] = 0;
110     fHistRhoVSEmbBkg[i] = 0;
111     fHistDeltaPtEmbArea[i] = 0;
112   }
113
114   SetMakeGeneralHistograms(kTRUE);
115 }
116
117 //________________________________________________________________________
118 AliAnalysisTaskDeltaPt::~AliAnalysisTaskDeltaPt()
119 {
120   // Destructor.
121 }
122
123 //________________________________________________________________________
124 void AliAnalysisTaskDeltaPt::UserCreateOutputObjects()
125 {
126   // Create user output.
127
128   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
129
130   fHistRCPhiEta = new TH2F("fHistRCPhiEta","Phi-Eta distribution of rigid cones", 50, -1, 1, 101, 0, TMath::Pi() * 2.02);
131   fHistRCPhiEta->GetXaxis()->SetTitle("#eta");
132   fHistRCPhiEta->GetYaxis()->SetTitle("#phi");
133   fOutput->Add(fHistRCPhiEta);
134
135   fHistRCPtExLJVSDPhiLJ = new TH2F("fHistRCPtExLJVSDPhiLJ","fHistRCPtExLJVSDPhiLJ", fNbins, fMinBinPt, fMaxBinPt, 128, -1.6, 4.8);
136   fHistRCPtExLJVSDPhiLJ->GetXaxis()->SetTitle("rigid cone #it{p}_{T} (GeV/#it{c})");
137   fHistRCPtExLJVSDPhiLJ->GetYaxis()->SetTitle("#Delta#phi");
138   fOutput->Add(fHistRCPtExLJVSDPhiLJ);
139
140   if (!fJetsName.IsNull()) {
141     fHistEmbJetPhiEta = new TH2F("fHistEmbJetPhiEta","Phi-Eta distribution of embedded jets", 50, -1, 1, 101, 0, TMath::Pi() * 2.02);
142     fHistEmbJetPhiEta->GetXaxis()->SetTitle("#eta");
143     fHistEmbJetPhiEta->GetYaxis()->SetTitle("#phi");
144     fOutput->Add(fHistEmbJetPhiEta);
145     
146     fHistEmbPartPhiEta = new TH2F("fHistEmbPartPhiEta","Phi-Eta distribution of embedded particles", 50, -1, 1, 101, 0, TMath::Pi() * 2.02);
147     fHistEmbPartPhiEta->GetXaxis()->SetTitle("#eta");
148     fHistEmbPartPhiEta->GetYaxis()->SetTitle("#phi");
149     fOutput->Add(fHistEmbPartPhiEta);
150   }
151
152   TString histname;
153
154   for (Int_t i = 0; i < 4; i++) {
155     histname = "fHistRCPt_";
156     histname += i;
157     fHistRCPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
158     fHistRCPt[i]->GetXaxis()->SetTitle("rigid cone #it{p}_{T} (GeV/#it{c})");
159     fHistRCPt[i]->GetYaxis()->SetTitle("counts");
160     fOutput->Add(fHistRCPt[i]);
161
162     histname = "fHistRCPtExLJ_";
163     histname += i;
164     fHistRCPtExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
165     fHistRCPtExLJ[i]->GetXaxis()->SetTitle("rigid cone #it{p}_{T}^{RC} (GeV/#it{c})");
166     fHistRCPtExLJ[i]->GetYaxis()->SetTitle("counts");
167     fOutput->Add(fHistRCPtExLJ[i]);
168
169     histname = "fHistRCPtRand_";
170     histname += i;
171     fHistRCPtRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
172     fHistRCPtRand[i]->GetXaxis()->SetTitle("rigid cone #it{p}_{T}^{RC} (GeV/#it{c})");
173     fHistRCPtRand[i]->GetYaxis()->SetTitle("counts");
174     fOutput->Add(fHistRCPtRand[i]);
175
176     histname = "fHistRhoVSRCPt_";
177     histname += i;
178     fHistRhoVSRCPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
179     fHistRhoVSRCPt[i]->GetXaxis()->SetTitle("A#rho (GeV/#it{c})");
180     fHistRhoVSRCPt[i]->GetYaxis()->SetTitle("rigid cone #it{p}_{T} (GeV/#it{c})");
181     fOutput->Add(fHistRhoVSRCPt[i]);
182
183     histname = "fHistDeltaPtRC_";
184     histname += i;
185     fHistDeltaPtRC[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
186     fHistDeltaPtRC[i]->GetXaxis()->SetTitle("#delta#it{p}_{T}^{RC} (GeV/#it{c})");
187     fHistDeltaPtRC[i]->GetYaxis()->SetTitle("counts");
188     fOutput->Add(fHistDeltaPtRC[i]);
189
190     histname = "fHistDeltaPtRCExLJ_";
191     histname += i;
192     fHistDeltaPtRCExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
193     fHistDeltaPtRCExLJ[i]->GetXaxis()->SetTitle("#delta#it{p}_{T}^{RC} (GeV/#it{c})");
194     fHistDeltaPtRCExLJ[i]->GetYaxis()->SetTitle("counts");
195     fOutput->Add(fHistDeltaPtRCExLJ[i]);
196
197     histname = "fHistDeltaPtRCRand_";
198     histname += i;
199     fHistDeltaPtRCRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
200     fHistDeltaPtRCRand[i]->GetXaxis()->SetTitle("#delta#it{p}_{T}^{RC} (GeV/#it{c})");
201     fHistDeltaPtRCRand[i]->GetYaxis()->SetTitle("counts");
202     fOutput->Add(fHistDeltaPtRCRand[i]);
203
204     if (!fEmbJetsName.IsNull()) {
205       histname = "fHistEmbJetsPtArea_";
206       histname += i;
207       fHistEmbJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
208       fHistEmbJetsPtArea[i]->GetXaxis()->SetTitle("area");
209       fHistEmbJetsPtArea[i]->GetYaxis()->SetTitle("embedded jet #it{p}_{T}^{raw} (GeV/#it{c})");
210       fOutput->Add(fHistEmbJetsPtArea[i]);
211
212       histname = "fHistEmbJetsCorrPtArea_";
213       histname += i;
214       fHistEmbJetsCorrPtArea[i] = new TH2F(histname.Data(), histname.Data(), 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins * 2, -fMaxBinPt, fMaxBinPt);
215       fHistEmbJetsCorrPtArea[i]->GetXaxis()->SetTitle("area");
216       fHistEmbJetsCorrPtArea[i]->GetYaxis()->SetTitle("embedded jet #it{p}_{T}^{corr} (GeV/#it{c})");
217       fOutput->Add(fHistEmbJetsCorrPtArea[i]);
218
219       histname = "fHistEmbPartPt_";
220       histname += i;
221       fHistEmbPartPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
222       fHistEmbPartPt[i]->GetXaxis()->SetTitle("embedded particle it{p}_{T}^{emb} (GeV/#it{c})");
223       fHistEmbPartPt[i]->GetYaxis()->SetTitle("counts");
224       fOutput->Add(fHistEmbPartPt[i]);
225
226       histname = "fHistDistEmbPartJetAxis_";
227       histname += i;
228       fHistDistEmbPartJetAxis[i] = new TH1F(histname.Data(), histname.Data(), 50, 0, 0.5);
229       fHistDistEmbPartJetAxis[i]->GetXaxis()->SetTitle("distance");
230       fHistDistEmbPartJetAxis[i]->GetYaxis()->SetTitle("counts");
231       fOutput->Add(fHistDistEmbPartJetAxis[i]);
232
233       histname = "fHistEmbNotFoundPhiEta_";
234       histname += i;
235       fHistEmbNotFoundPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 40, -2, 2, 64, 0, 6.4);
236       fHistEmbNotFoundPhiEta[i]->GetXaxis()->SetTitle("#eta");
237       fHistEmbNotFoundPhiEta[i]->GetYaxis()->SetTitle("#phi");
238       fOutput->Add(fHistEmbNotFoundPhiEta[i]);
239
240       histname = "fHistEmbBkgArea_";
241       histname += i;
242       fHistEmbBkgArea[i] = new TH2F(histname.Data(), histname.Data(), 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
243       fHistEmbBkgArea[i]->GetXaxis()->SetTitle("area");
244       fHistEmbBkgArea[i]->GetYaxis()->SetTitle("background of embedded track (GeV/#it{c})");
245       fOutput->Add(fHistEmbBkgArea[i]);
246
247       histname = "fHistRhoVSEmbBkg_";
248       histname += i;
249       fHistRhoVSEmbBkg[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
250       fHistRhoVSEmbBkg[i]->GetXaxis()->SetTitle("A#rho (GeV/#it{c})");
251       fHistRhoVSEmbBkg[i]->GetYaxis()->SetTitle("background of embedded track (GeV/#it{c})");
252       fOutput->Add(fHistRhoVSEmbBkg[i]);
253       
254       histname = "fHistDeltaPtEmbArea_";
255       histname += i;
256       fHistDeltaPtEmbArea[i] = new TH2F(histname.Data(), histname.Data(), 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins * 2, -fMaxBinPt, fMaxBinPt);
257       fHistDeltaPtEmbArea[i]->GetXaxis()->SetTitle("area");
258       fHistDeltaPtEmbArea[i]->GetYaxis()->SetTitle("#delta#it{p}_{T}^{emb} (GeV/#it{c})");
259       fOutput->Add(fHistDeltaPtEmbArea[i]);
260     }
261   }
262
263   PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
264 }
265
266 //________________________________________________________________________
267 Bool_t AliAnalysisTaskDeltaPt::FillHistograms()
268 {
269   // Fill histograms.
270
271   Int_t *sortedJets = GetSortedArray(fJets);
272   
273   AliEmcalJet* jet = 0;
274
275   if (sortedJets && sortedJets[0] > 0) 
276     jet = static_cast<AliEmcalJet*>(fJets->At(sortedJets[0]));
277
278   // ************
279   // Random cones
280   // _________________________________
281   
282   const Float_t rcArea = fJetRadius * fJetRadius * TMath::Pi();
283
284   for (Int_t i = 0; i < fRCperEvent; i++) {
285     // Simple random cones
286     Float_t RCpt = 0;
287     Float_t RCeta = 0;
288     Float_t RCphi = 0;
289     GetRandomCone(RCpt, RCeta, RCphi, 0);
290     if (RCpt > 0) {
291       fHistRCPhiEta->Fill(RCeta, RCphi);
292       fHistRhoVSRCPt[fCentBin]->Fill(fRhoVal * rcArea, RCpt);
293
294       fHistRCPt[fCentBin]->Fill(RCpt);
295       fHistDeltaPtRC[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
296     }
297   
298     // Random cones far from leading jet
299     RCpt = 0;
300     RCeta = 0;
301     RCphi = 0;
302     GetRandomCone(RCpt, RCeta, RCphi, jet);
303     if (RCpt > 0) {
304       if (jet) {
305         Float_t dphi = RCphi - jet->Phi();
306         if (dphi > 4.8) dphi -= TMath::Pi() * 2;
307         if (dphi < -1.6) dphi += TMath::Pi() * 2; 
308         fHistRCPtExLJVSDPhiLJ->Fill(RCpt, dphi);
309       }
310       fHistRCPtExLJ[fCentBin]->Fill(RCpt);
311       fHistDeltaPtRCExLJ[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
312     }
313     
314     // Random cones with randomized particles
315     RCpt = 0;
316     RCeta = 0;
317     RCphi = 0;
318     GetRandomCone(RCpt, RCeta, RCphi, 0, fRandTracks, fRandCaloClusters);
319     if (RCpt > 0) {
320       fHistRCPtRand[fCentBin]->Fill(RCpt);
321       fHistDeltaPtRCRand[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
322     }  
323   }
324   
325   // ************
326   // Embedding
327   // _________________________________
328
329   if (!fEmbJets)
330     return kTRUE;
331
332   AliEmcalJet *embJet  = 0;
333   TObject     *embPart = 0;
334
335   DoEmbJetLoop(embJet, embPart);
336
337   if (embJet) {
338     Double_t probePt = 0;
339     Double_t probeEta = 0;
340     Double_t probePhi = 0;
341
342     AliVCluster *cluster = dynamic_cast<AliVCluster*>(embPart);
343     if (cluster) {
344       TLorentzVector nPart;
345       cluster->GetMomentum(nPart, fVertex);
346
347       probeEta = nPart.Eta();
348       probePhi = nPart.Phi();
349       probePt = nPart.Pt();
350     }
351     else {
352       AliVTrack *track = dynamic_cast<AliVTrack*>(embPart);
353       if (track) {
354         probeEta = track->Eta();
355         probePhi = track->Phi();
356         probePt = track->Pt();
357       }
358       else {
359         AliWarning(Form("%s - Embedded jet found but embedded particle not found (wrong type?)!", GetName()));
360         return kTRUE;
361       }
362     }
363     Double_t distProbeJet = TMath::Sqrt((embJet->Eta() - probeEta) * (embJet->Eta() - probeEta) + (embJet->Phi() - probePhi) * (embJet->Phi() - probePhi));
364
365     fHistEmbPartPt[fCentBin]->Fill(probePt);
366     fHistEmbPartPhiEta->Fill(probeEta, probePhi);
367     
368     fHistDistEmbPartJetAxis[fCentBin]->Fill(distProbeJet);
369
370     fHistEmbJetsPtArea[fCentBin]->Fill(embJet->Area(), embJet->Pt());
371     fHistEmbJetsCorrPtArea[fCentBin]->Fill(embJet->Area(), embJet->Pt() - fRhoVal * embJet->Area());
372     fHistEmbJetPhiEta->Fill(embJet->Eta(), embJet->Phi());
373
374     fHistEmbBkgArea[fCentBin]->Fill(embJet->Area(), embJet->Pt() - probePt);
375     fHistRhoVSEmbBkg[fCentBin]->Fill(fRhoVal * embJet->Area(), embJet->Pt() - probePt);
376     fHistDeltaPtEmbArea[fCentBin]->Fill(embJet->Area(), embJet->Pt() - embJet->Area() * fRhoVal - probePt);
377
378   }
379   else {
380     if (fEmbTracks)
381       DoEmbTrackLoop();
382     if (fEmbCaloClusters)
383       DoEmbClusterLoop();
384     if (fEmbeddedTrackId >= 0) {
385       AliVTrack *track2 = static_cast<AliVTrack*>(fEmbTracks->At(fEmbeddedTrackId));
386       fHistEmbNotFoundPhiEta[fCentBin]->Fill(track2->Eta(), track2->Phi());
387     }
388     else if (fEmbeddedClusterId >= 0) {
389       AliVCluster *cluster2 = static_cast<AliVCluster*>(fEmbCaloClusters->At(fEmbeddedClusterId));
390       TLorentzVector nPart;
391       cluster2->GetMomentum(nPart, fVertex);
392       fHistEmbNotFoundPhiEta[fCentBin]->Fill(nPart.Eta(), nPart.Phi());
393     }
394     else {
395       AliWarning(Form("%s - Embedded particle not found!", GetName()));
396     }
397   }
398
399   return kTRUE;
400 }
401
402 //________________________________________________________________________
403 void AliAnalysisTaskDeltaPt::DoEmbTrackLoop()
404 {
405   // Do track loop.
406
407   if (!fEmbTracks)
408     return;
409
410   fEmbeddedTrackId = -1;
411
412   Int_t ntracks = fEmbTracks->GetEntriesFast();
413
414   for (Int_t i = 0; i < ntracks; i++) {
415
416     AliVParticle* track = static_cast<AliVParticle*>(fEmbTracks->At(i)); // pointer to reconstructed to track  
417
418     if (!track) {
419       AliError(Form("Could not retrieve track %d",i)); 
420       continue; 
421     }
422
423     AliVTrack* vtrack = dynamic_cast<AliVTrack*>(track); 
424     
425     if (vtrack && !AcceptTrack(vtrack, kTRUE)) 
426       continue;
427
428     if (track->GetLabel() == 100) {
429       fEmbeddedTrackId = i;
430       continue;
431     }
432   }
433 }
434
435 //________________________________________________________________________
436 void AliAnalysisTaskDeltaPt::DoEmbClusterLoop()
437 {
438   // Do cluster loop.
439
440   if (!fEmbCaloClusters)
441     return;
442
443   fEmbeddedClusterId = -1;
444
445   Int_t nclusters =  fEmbCaloClusters->GetEntriesFast();
446
447   for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
448     AliVCluster* cluster = static_cast<AliVCluster*>(fEmbCaloClusters->At(iClusters));
449     if (!cluster) {
450       AliError(Form("Could not receive cluster %d", iClusters));
451       continue;
452     }  
453
454     if (!AcceptCluster(cluster, kTRUE)) 
455       continue;
456
457     if (cluster->Chi2() == 100) {
458       fEmbeddedClusterId = iClusters;
459       continue;
460     }
461   }
462 }
463
464 //________________________________________________________________________
465 void AliAnalysisTaskDeltaPt::DoEmbJetLoop(AliEmcalJet* &embJet, TObject* &embPart)
466 {
467   // Do the embedded jet loop.
468
469   if (!fEmbJets)
470     return;
471
472   TLorentzVector maxClusVect;
473
474   Int_t nembjets = fEmbJets->GetEntriesFast();
475
476   for (Int_t ij = 0; ij < nembjets; ij++) {
477       
478     AliEmcalJet* jet = static_cast<AliEmcalJet*>(fEmbJets->At(ij));
479       
480     if (!jet) {
481       AliError(Form("Could not receive jet %d", ij));
482       continue;
483     } 
484       
485     if (!AcceptJet(jet))
486       continue;
487
488     if (!jet->IsMC())
489       continue;
490
491     AliVParticle *maxTrack = 0;
492
493     for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
494       AliVParticle *track = jet->TrackAt(it, fEmbTracks);
495       
496       if (!track) 
497         continue;
498
499       if (track->GetLabel() != 100)
500         continue;
501       
502       if (!maxTrack || track->Pt() > maxTrack->Pt())
503         maxTrack = track;
504     }
505     
506     AliVCluster *maxClus = 0;
507
508     for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
509       AliVCluster *cluster = jet->ClusterAt(ic, fEmbCaloClusters);
510       
511       if (!cluster) 
512         continue;
513
514       if (cluster->Chi2() != 100)
515         continue;
516       
517       TLorentzVector nPart;
518       cluster->GetMomentum(nPart, fVertex);
519       
520       if (!maxClus || nPart.Et() > maxClusVect.Et()) {
521         new (&maxClusVect) TLorentzVector(nPart);
522         maxClus = cluster;
523       } 
524     }
525
526     if ((maxClus && maxTrack && maxClusVect.Et() > maxTrack->Pt()) || (maxClus && !maxTrack)) {
527       embPart = maxClus;
528       embJet = jet;
529     }
530     else if (maxTrack) {
531       embPart = maxTrack;
532       embJet = jet;
533     }
534
535     return;  // MC jet found, exit
536   }
537 }
538
539 //________________________________________________________________________
540 void AliAnalysisTaskDeltaPt::GetRandomCone(Float_t &pt, Float_t &eta, Float_t &phi,
541                                         AliEmcalJet *jet, TClonesArray* tracks, TClonesArray* clusters) const
542 {
543   // Get rigid cone.
544
545   if (!tracks)
546     tracks = fTracks;
547
548   if (!clusters)
549     clusters = fCaloClusters;
550
551   eta = 0;
552   phi = 0;
553   pt = 0;
554
555   if (!tracks && !clusters)
556     return;
557
558   Float_t LJeta = 999;
559   Float_t LJphi = 999;
560
561   if (jet) {
562     LJeta = jet->Eta();
563     LJphi = jet->Phi();
564   }
565
566   Float_t maxEta = fJetMaxEta;
567   Float_t minEta = fJetMinEta;
568   Float_t maxPhi = fJetMaxPhi;
569   Float_t minPhi = fJetMinPhi;
570
571   if (maxPhi > TMath::Pi() * 2) maxPhi = TMath::Pi() * 2;
572   if (minPhi < 0) minPhi = 0;
573   
574   Float_t dLJ = 0;
575   Int_t repeats = 0;
576   do {
577     eta = gRandom->Rndm() * (maxEta - minEta) + minEta;
578     phi = gRandom->Rndm() * (maxPhi - minPhi) + minPhi;
579     dLJ = TMath::Sqrt((LJeta - eta) * (LJeta - eta) + (LJphi - phi) * (LJphi - phi));
580     repeats++;
581   } while (dLJ < fMinRC2LJ && repeats < 999);
582
583   if (repeats == 999) {
584     AliWarning(Form("%s: Could not get random cone!", GetName()));
585     return;
586   }
587
588   if (clusters) {
589     Int_t nclusters =  clusters->GetEntriesFast();
590     for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
591       AliVCluster* cluster = static_cast<AliVCluster*>(clusters->At(iClusters));
592       if (!cluster) {
593         AliError(Form("Could not receive cluster %d", iClusters));
594         continue;
595       }  
596       
597       if (!AcceptCluster(cluster, fMCAna))
598         continue;
599       
600       TLorentzVector nPart;
601       cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
602      
603       Float_t d = TMath::Sqrt((nPart.Eta() - eta) * (nPart.Eta() - eta) + (nPart.Phi() - phi) * (nPart.Phi() - phi));
604
605       if (d <= fJetRadius) 
606         pt += nPart.Pt();
607     }
608   }
609
610   if (tracks) {
611     Int_t ntracks = tracks->GetEntriesFast();
612     for(Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
613       AliVTrack* track = static_cast<AliVTrack*>(tracks->At(iTracks));         
614       if(!track) {
615         AliError(Form("Could not retrieve track %d",iTracks)); 
616         continue; 
617       }
618       
619       if (!AcceptTrack(track, fMCAna)) 
620         continue;
621       
622       Float_t tracketa = track->Eta();
623       Float_t trackphi = track->Phi();
624       
625       if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi + 2 * TMath::Pi()))
626         trackphi += 2 * TMath::Pi();
627       if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi - 2 * TMath::Pi()))
628         trackphi -= 2 * TMath::Pi();
629       
630       Float_t d = TMath::Sqrt((tracketa - eta) * (tracketa - eta) + (trackphi - phi) * (trackphi - phi));
631       if (d <= fJetRadius)
632         pt += track->Pt();
633     }
634   }
635 }
636
637 //________________________________________________________________________
638 void AliAnalysisTaskDeltaPt::ExecOnce()
639 {
640   // Initialize the analysis.
641
642   if (!fEmbJetsName.IsNull() && !fEmbJets) {
643     fEmbJets =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbJetsName));
644     if (!fEmbJets) {
645       AliError(Form("%s: Could not retrieve embedded jets %s!", GetName(), fEmbJetsName.Data()));
646       return;
647     }
648     else if (!fEmbJets->GetClass()->GetBaseClass("AliEmcalJet")) {
649       AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fEmbJetsName.Data())); 
650       fEmbJets = 0;
651       return;
652     }
653   }
654
655   if (!fEmbCaloName.IsNull() && !fEmbCaloClusters) {
656     fEmbCaloClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbCaloName));
657     if (!fEmbCaloClusters) {
658       AliError(Form("%s: Could not retrieve embedded clusters %s!", GetName(), fEmbCaloName.Data()));
659       return;
660     }
661     else if (!fEmbCaloClusters->GetClass()->GetBaseClass("AliVCluster") && !fEmbCaloClusters->GetClass()->GetBaseClass("AliEmcalParticle")) {
662       AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fEmbCaloName.Data())); 
663       fEmbCaloClusters = 0;
664       return;
665     }
666   }
667
668   if (!fEmbTracksName.IsNull() && !fEmbTracks) {
669     fEmbTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbTracksName));
670     if (!fEmbTracks) {
671       AliError(Form("%s: Could not retrieve embedded tracks %s!", GetName(), fEmbTracksName.Data()));
672       return;
673     }
674     else if (!fEmbTracks->GetClass()->GetBaseClass("AliVParticle") && !fEmbTracks->GetClass()->GetBaseClass("AliEmcalParticle")) {
675       AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fEmbTracksName.Data())); 
676       fEmbTracks = 0;
677       return;
678     }
679   }
680
681   if (!fRandCaloName.IsNull() && !fRandCaloClusters) {
682     fRandCaloClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fRandCaloName));
683     if (!fRandCaloClusters) {
684       AliError(Form("%s: Could not retrieve randomized clusters %s!", GetName(), fRandCaloName.Data()));
685       return;
686     }
687     else if (!fRandCaloClusters->GetClass()->GetBaseClass("AliVCluster") && !fRandCaloClusters->GetClass()->GetBaseClass("AliEmcalParticle")) {
688       AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fRandCaloName.Data())); 
689       fRandCaloClusters = 0;
690       return;
691     }
692   }
693
694   if (!fRandTracksName.IsNull() && !fRandTracks) {
695     fRandTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fRandTracksName));
696     if (!fRandTracks) {
697       AliError(Form("%s: Could not retrieve randomized tracks %s!", GetName(), fRandTracksName.Data()));
698       return;
699     }
700     else if (!fRandTracks->GetClass()->GetBaseClass("AliVParticle") && !fRandTracks->GetClass()->GetBaseClass("AliEmcalParticle")) {
701       AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fRandTracksName.Data())); 
702       fRandTracks = 0;
703       return;
704     }
705   }
706
707   AliAnalysisTaskEmcalJet::ExecOnce();
708
709   if (fRCperEvent < 0) {
710     Double_t area = (fJetMaxEta - fJetMinEta) * (fJetMaxPhi - fJetMinPhi);
711     Double_t jetArea = TMath::Pi() * fJetRadius * fJetRadius;
712     fRCperEvent = TMath::FloorNint(area / jetArea - 0.5);
713     if (fRCperEvent == 0)
714       fRCperEvent = 1;
715   }
716
717   if (fMinRC2LJ < 0)
718     fMinRC2LJ = fJetRadius * 1.5;
719
720   const Float_t maxDist = TMath::Max(fJetMaxPhi - fJetMinPhi, fJetMaxEta - fJetMinEta) / 2;
721   if (fMinRC2LJ > maxDist) {
722     AliWarning(Form("The parameter fMinRC2LJ = %f is too large for the considered acceptance. "
723                     "Will use fMinRC2LJ = %f", fMinRC2LJ, maxDist));
724     fMinRC2LJ = maxDist;
725   }
726 }
727
728 //________________________________________________________________________
729 void AliAnalysisTaskDeltaPt::Terminate(Option_t *) 
730 {
731   // Called once at the end of the analysis.
732 }