8f9a6d372d54a3b955c2a1cd0b672f28c5192fd1
[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     fHistEmbJetsPt[i] = 0;
60     fHistEmbJetsCorrPt[i] = 0;
61     fHistEmbJetsArea[i] = 0;
62     fHistEmbPartPt[i] = 0;
63     fHistDistEmbPartJetAxis[i] = 0;
64     fHistRhoVSEmbBkg[i] = 0;
65     fHistDeltaPtEmb[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     fHistEmbJetsPt[i] = 0;
106     fHistEmbJetsCorrPt[i] = 0;
107     fHistEmbJetsArea[i] = 0;
108     fHistEmbPartPt[i] = 0;
109     fHistDistEmbPartJetAxis[i] = 0;
110     fHistRhoVSEmbBkg[i] = 0;
111     fHistDeltaPtEmb[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 p_{T} [GeV/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 p_{T} [GeV/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 p_{T}^{RC} [GeV/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 p_{T}^{RC} [GeV/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*2, fNbins, fMinBinPt, fMaxBinPt);
179     fHistRhoVSRCPt[i]->GetXaxis()->SetTitle("A#rho [GeV/c]");
180     fHistRhoVSRCPt[i]->GetYaxis()->SetTitle("rigid cone p_{T} [GeV/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("#deltap_{T}^{RC} [GeV/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("#deltap_{T}^{RC} [GeV/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("#deltap_{T}^{RC} [GeV/c]");
201     fHistDeltaPtRCRand[i]->GetYaxis()->SetTitle("counts");
202     fOutput->Add(fHistDeltaPtRCRand[i]);
203
204     if (!fEmbJetsName.IsNull()) {
205       histname = "fHistEmbJetsPt_";
206       histname += i;
207       fHistEmbJetsPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
208       fHistEmbJetsPt[i]->GetXaxis()->SetTitle("embedded jet p_{T}^{raw} [GeV/c]");
209       fHistEmbJetsPt[i]->GetYaxis()->SetTitle("counts");
210       fOutput->Add(fHistEmbJetsPt[i]);
211
212       histname = "fHistEmbJetsCorrPt_";
213       histname += i;
214       fHistEmbJetsCorrPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
215       fHistEmbJetsCorrPt[i]->GetXaxis()->SetTitle("embedded jet p_{T}^{corr} [GeV/c]");
216       fHistEmbJetsCorrPt[i]->GetYaxis()->SetTitle("counts");
217       fOutput->Add(fHistEmbJetsCorrPt[i]);
218
219       histname = "fHistEmbJetsArea_";
220       histname += i;
221       fHistEmbJetsArea[i] = new TH1F(histname.Data(), histname.Data(), 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3);
222       fHistEmbJetsArea[i]->GetXaxis()->SetTitle("area");
223       fHistEmbJetsArea[i]->GetYaxis()->SetTitle("counts");
224       fOutput->Add(fHistEmbJetsArea[i]);
225
226       histname = "fHistEmbPartPt_";
227       histname += i;
228       fHistEmbPartPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
229       fHistEmbPartPt[i]->GetXaxis()->SetTitle("embedded particle p_{T}^{emb} [GeV/c]");
230       fHistEmbPartPt[i]->GetYaxis()->SetTitle("counts");
231       fOutput->Add(fHistEmbPartPt[i]);
232
233       histname = "fHistDistEmbPartJetAxis_";
234       histname += i;
235       fHistDistEmbPartJetAxis[i] = new TH1F(histname.Data(), histname.Data(), 50, 0, 0.5);
236       fHistDistEmbPartJetAxis[i]->GetXaxis()->SetTitle("distance");
237       fHistDistEmbPartJetAxis[i]->GetYaxis()->SetTitle("counts");
238       fOutput->Add(fHistDistEmbPartJetAxis[i]);
239
240       histname = "fHistEmbNotFoundPhiEta_";
241       histname += i;
242       fHistEmbNotFoundPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 40, -2, 2, 64, 0, 6.4);
243       fHistEmbNotFoundPhiEta[i]->GetXaxis()->SetTitle("#eta");
244       fHistEmbNotFoundPhiEta[i]->GetYaxis()->SetTitle("#phi");
245       fOutput->Add(fHistEmbNotFoundPhiEta[i]);
246
247       histname = "fHistRhoVSEmbBkg_";
248       histname += i;
249       fHistRhoVSEmbBkg[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt*2, fNbins, fMinBinPt, fMaxBinPt);
250       fHistRhoVSEmbBkg[i]->GetXaxis()->SetTitle("A#rho [GeV/c]");
251       fHistRhoVSEmbBkg[i]->GetYaxis()->SetTitle("background of embedded track [GeV/c]");
252       fOutput->Add(fHistRhoVSEmbBkg[i]);
253       
254       histname = "fHistDeltaPtEmb_";
255       histname += i;
256       fHistDeltaPtEmb[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
257       fHistDeltaPtEmb[i]->GetXaxis()->SetTitle("#deltap_{T}^{emb} [GeV/c]");
258       fHistDeltaPtEmb[i]->GetYaxis()->SetTitle("counts");
259       fOutput->Add(fHistDeltaPtEmb[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   // Simple random cones
285   for (Int_t i = 0; i < fRCperEvent; i++) {
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, RCpt);
293
294       fHistRCPt[fCentBin]->Fill(RCpt / rcArea);
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 / rcArea);
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 / rcArea);
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     fHistEmbJetsPt[fCentBin]->Fill(embJet->Pt());
371     fHistEmbJetsCorrPt[fCentBin]->Fill(embJet->Pt() - fRhoVal * embJet->Area());
372     fHistEmbJetsArea[fCentBin]->Fill(embJet->Area());
373     fHistEmbJetPhiEta->Fill(embJet->Eta(), embJet->Phi());
374
375     fHistDeltaPtEmb[fCentBin]->Fill(embJet->Pt() - embJet->Area() * fRhoVal - probePt);
376     fHistRhoVSEmbBkg[fCentBin]->Fill(fRhoVal, embJet->Pt() - probePt);
377   }
378   else {
379     if (fEmbTracks)
380       DoEmbTrackLoop();
381     if (fEmbCaloClusters)
382       DoEmbClusterLoop();
383     if (fEmbeddedTrackId >= 0) {
384       AliVTrack *track2 = static_cast<AliVTrack*>(fEmbTracks->At(fEmbeddedTrackId));
385       fHistEmbNotFoundPhiEta[fCentBin]->Fill(track2->Eta(), track2->Phi());
386     }
387     else if (fEmbeddedClusterId >= 0) {
388       AliVCluster *cluster2 = static_cast<AliVCluster*>(fEmbCaloClusters->At(fEmbeddedClusterId));
389       TLorentzVector nPart;
390       cluster2->GetMomentum(nPart, fVertex);
391       fHistEmbNotFoundPhiEta[fCentBin]->Fill(nPart.Eta(), nPart.Phi());
392     }
393     else {
394       AliWarning(Form("%s - Embedded particle not found!", GetName()));
395     }
396   }
397
398   return kTRUE;
399 }
400
401 //________________________________________________________________________
402 void AliAnalysisTaskDeltaPt::DoEmbTrackLoop()
403 {
404   // Do track loop.
405
406   if (!fEmbTracks)
407     return;
408
409   fEmbeddedTrackId = -1;
410
411   Int_t ntracks = fEmbTracks->GetEntriesFast();
412
413   for (Int_t i = 0; i < ntracks; i++) {
414
415     AliVParticle* track = static_cast<AliVParticle*>(fEmbTracks->At(i)); // pointer to reconstructed to track  
416
417     if (!track) {
418       AliError(Form("Could not retrieve track %d",i)); 
419       continue; 
420     }
421
422     AliVTrack* vtrack = dynamic_cast<AliVTrack*>(track); 
423     
424     if (vtrack && !AcceptTrack(vtrack, kTRUE)) 
425       continue;
426
427     if (track->GetLabel() == 100) {
428       fEmbeddedTrackId = i;
429       continue;
430     }
431   }
432 }
433
434 //________________________________________________________________________
435 void AliAnalysisTaskDeltaPt::DoEmbClusterLoop()
436 {
437   // Do cluster loop.
438
439   if (!fEmbCaloClusters)
440     return;
441
442   fEmbeddedClusterId = -1;
443
444   Int_t nclusters =  fEmbCaloClusters->GetEntriesFast();
445
446   for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
447     AliVCluster* cluster = static_cast<AliVCluster*>(fEmbCaloClusters->At(iClusters));
448     if (!cluster) {
449       AliError(Form("Could not receive cluster %d", iClusters));
450       continue;
451     }  
452
453     if (!AcceptCluster(cluster, kTRUE)) 
454       continue;
455
456     if (cluster->Chi2() == 100) {
457       fEmbeddedClusterId = iClusters;
458       continue;
459     }
460   }
461 }
462
463 //________________________________________________________________________
464 void AliAnalysisTaskDeltaPt::DoEmbJetLoop(AliEmcalJet* &embJet, TObject* &embPart)
465 {
466   // Do the embedded jet loop.
467
468   if (!fEmbJets)
469     return;
470
471   TLorentzVector maxClusVect;
472
473   Int_t nembjets = fEmbJets->GetEntriesFast();
474
475   for (Int_t ij = 0; ij < nembjets; ij++) {
476       
477     AliEmcalJet* jet = static_cast<AliEmcalJet*>(fEmbJets->At(ij));
478       
479     if (!jet) {
480       AliError(Form("Could not receive jet %d", ij));
481       continue;
482     } 
483       
484     if (!AcceptJet(jet))
485       continue;
486
487     if (!jet->IsMC())
488       continue;
489
490     AliVParticle *maxTrack = 0;
491
492     for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
493       AliVParticle *track = jet->TrackAt(it, fEmbTracks);
494       
495       if (!track) 
496         continue;
497
498       if (track->GetLabel() != 100)
499         continue;
500       
501       if (!maxTrack || track->Pt() > maxTrack->Pt())
502         maxTrack = track;
503     }
504     
505     AliVCluster *maxClus = 0;
506
507     for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
508       AliVCluster *cluster = jet->ClusterAt(ic, fEmbCaloClusters);
509       
510       if (!cluster) 
511         continue;
512
513       if (cluster->Chi2() != 100)
514         continue;
515       
516       TLorentzVector nPart;
517       cluster->GetMomentum(nPart, fVertex);
518       
519       if (!maxClus || nPart.Et() > maxClusVect.Et()) {
520         new (&maxClusVect) TLorentzVector(nPart);
521         maxClus = cluster;
522       } 
523     }
524
525     if ((maxClus && maxTrack && maxClusVect.Et() > maxTrack->Pt()) || (maxClus && !maxTrack)) {
526       embPart = maxClus;
527       embJet = jet;
528     }
529     else if (maxTrack) {
530       embPart = maxTrack;
531       embJet = jet;
532     }
533
534     return;  // MC jet found, exit
535   }
536 }
537
538 //________________________________________________________________________
539 void AliAnalysisTaskDeltaPt::GetRandomCone(Float_t &pt, Float_t &eta, Float_t &phi,
540                                         AliEmcalJet *jet, TClonesArray* tracks, TClonesArray* clusters) const
541 {
542   // Get rigid cone.
543
544   if (!tracks)
545     tracks = fTracks;
546
547   if (!clusters)
548     clusters = fCaloClusters;
549
550   eta = 0;
551   phi = 0;
552   pt = 0;
553
554   if (!tracks && !clusters)
555     return;
556
557   Float_t LJeta = 999;
558   Float_t LJphi = 999;
559
560   if (jet) {
561     LJeta = jet->Eta();
562     LJphi = jet->Phi();
563   }
564
565   Float_t maxEta = fJetMaxEta;
566   Float_t minEta = fJetMinEta;
567   Float_t maxPhi = fJetMaxPhi;
568   Float_t minPhi = fJetMinPhi;
569
570   if (maxPhi > TMath::Pi() * 2) maxPhi = TMath::Pi() * 2;
571   if (minPhi < 0) minPhi = 0;
572   
573   Float_t dLJ = 0;
574   Int_t repeats = 0;
575   do {
576     eta = gRandom->Rndm() * (maxEta - minEta) + minEta;
577     phi = gRandom->Rndm() * (maxPhi - minPhi) + minPhi;
578     dLJ = TMath::Sqrt((LJeta - eta) * (LJeta - eta) + (LJphi - phi) * (LJphi - phi));
579     repeats++;
580   } while (dLJ < fMinRC2LJ && repeats < 999);
581
582   if (repeats == 999) {
583     AliWarning(Form("%s: Could not get random cone!", GetName()));
584     return;
585   }
586
587   if (clusters) {
588     Int_t nclusters =  clusters->GetEntriesFast();
589     for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
590       AliVCluster* cluster = static_cast<AliVCluster*>(clusters->At(iClusters));
591       if (!cluster) {
592         AliError(Form("Could not receive cluster %d", iClusters));
593         continue;
594       }  
595       
596       if (!AcceptCluster(cluster, fMCAna))
597         continue;
598       
599       TLorentzVector nPart;
600       cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
601      
602       Float_t d = TMath::Sqrt((nPart.Eta() - eta) * (nPart.Eta() - eta) + (nPart.Phi() - phi) * (nPart.Phi() - phi));
603
604       if (d <= fJetRadius) 
605         pt += nPart.Pt();
606     }
607   }
608
609   if (tracks) {
610     Int_t ntracks = tracks->GetEntriesFast();
611     for(Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
612       AliVTrack* track = static_cast<AliVTrack*>(tracks->At(iTracks));         
613       if(!track) {
614         AliError(Form("Could not retrieve track %d",iTracks)); 
615         continue; 
616       }
617       
618       if (!AcceptTrack(track, fMCAna)) 
619         continue;
620       
621       Float_t tracketa = track->Eta();
622       Float_t trackphi = track->Phi();
623       
624       if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi + 2 * TMath::Pi()))
625         trackphi += 2 * TMath::Pi();
626       if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi - 2 * TMath::Pi()))
627         trackphi -= 2 * TMath::Pi();
628       
629       Float_t d = TMath::Sqrt((tracketa - eta) * (tracketa - eta) + (trackphi - phi) * (trackphi - phi));
630       if (d <= fJetRadius)
631         pt += track->Pt();
632     }
633   }
634 }
635
636 //________________________________________________________________________
637 void AliAnalysisTaskDeltaPt::ExecOnce()
638 {
639   // Initialize the analysis.
640
641   if (!fEmbJetsName.IsNull() && !fEmbJets) {
642     fEmbJets =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbJetsName));
643     if (!fEmbJets) {
644       AliError(Form("%s: Could not retrieve embedded jets %s!", GetName(), fEmbJetsName.Data()));
645       return;
646     }
647     else if (!fEmbJets->GetClass()->GetBaseClass("AliEmcalJet")) {
648       AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fEmbJetsName.Data())); 
649       fEmbJets = 0;
650       return;
651     }
652   }
653
654   if (!fEmbCaloName.IsNull() && !fEmbCaloClusters) {
655     fEmbCaloClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbCaloName));
656     if (!fEmbCaloClusters) {
657       AliError(Form("%s: Could not retrieve embedded clusters %s!", GetName(), fEmbCaloName.Data()));
658       return;
659     }
660     else if (!fEmbCaloClusters->GetClass()->GetBaseClass("AliVCluster") && !fEmbCaloClusters->GetClass()->GetBaseClass("AliEmcalParticle")) {
661       AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fEmbCaloName.Data())); 
662       fEmbCaloClusters = 0;
663       return;
664     }
665   }
666
667   if (!fEmbTracksName.IsNull() && !fEmbTracks) {
668     fEmbTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbTracksName));
669     if (!fEmbTracks) {
670       AliError(Form("%s: Could not retrieve embedded tracks %s!", GetName(), fEmbTracksName.Data()));
671       return;
672     }
673     else if (!fEmbTracks->GetClass()->GetBaseClass("AliVParticle") && !fEmbTracks->GetClass()->GetBaseClass("AliEmcalParticle")) {
674       AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fEmbTracksName.Data())); 
675       fEmbTracks = 0;
676       return;
677     }
678   }
679
680   if (!fRandCaloName.IsNull() && !fRandCaloClusters) {
681     fRandCaloClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fRandCaloName));
682     if (!fRandCaloClusters) {
683       AliError(Form("%s: Could not retrieve randomized clusters %s!", GetName(), fRandCaloName.Data()));
684       return;
685     }
686     else if (!fRandCaloClusters->GetClass()->GetBaseClass("AliVCluster") && !fRandCaloClusters->GetClass()->GetBaseClass("AliEmcalParticle")) {
687       AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fRandCaloName.Data())); 
688       fRandCaloClusters = 0;
689       return;
690     }
691   }
692
693   if (!fRandTracksName.IsNull() && !fRandTracks) {
694     fRandTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fRandTracksName));
695     if (!fRandTracks) {
696       AliError(Form("%s: Could not retrieve randomized tracks %s!", GetName(), fRandTracksName.Data()));
697       return;
698     }
699     else if (!fRandTracks->GetClass()->GetBaseClass("AliVParticle") && !fRandTracks->GetClass()->GetBaseClass("AliEmcalParticle")) {
700       AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fRandTracksName.Data())); 
701       fRandTracks = 0;
702       return;
703     }
704   }
705
706   AliAnalysisTaskEmcalJet::ExecOnce();
707
708   if (fRCperEvent < 0) {
709     Double_t area = (fJetMaxEta - fJetMinEta) * (fJetMaxPhi - fJetMinPhi);
710     Double_t jetArea = TMath::Pi() * fJetRadius * fJetRadius;
711     fRCperEvent = TMath::FloorNint(area / jetArea - 0.5);
712     if (fRCperEvent == 0)
713       fRCperEvent = 1;
714   }
715
716   if (fMinRC2LJ < 0)
717     fMinRC2LJ = fJetRadius * 1.5;
718
719   const Float_t maxDist = TMath::Max(fJetMaxPhi - fJetMinPhi, fJetMaxEta - fJetMinEta) / 2;
720   if (fMinRC2LJ > maxDist) {
721     AliWarning(Form("The parameter fMinRC2LJ = %f is too large for the considered acceptance. "
722                     "Will use fMinRC2LJ = %f", fMinRC2LJ, maxDist));
723     fMinRC2LJ = maxDist;
724   }
725 }
726
727 //________________________________________________________________________
728 void AliAnalysisTaskDeltaPt::Terminate(Option_t *) 
729 {
730   // Called once at the end of the analysis.
731 }