282bc7740f477b0595856236d96e35c2ef55de67
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskEmcalDiJetBase.cxx
1 //
2 // Dijet base analysis task.
3 //
4 // Author: M.Verweij
5
6 #include <TClonesArray.h>
7 #include <TH1F.h>
8 #include <TH2F.h>
9 #include <TH3F.h>
10 #include <THnSparse.h>
11 #include <TList.h>
12 #include <TLorentzVector.h>
13 #include <TProfile.h>
14 #include <TChain.h>
15 #include <TSystem.h>
16 #include <TFile.h>
17 #include <TKey.h>
18
19 #include "AliVCluster.h"
20 #include "AliVTrack.h"
21 #include "AliEmcalJet.h"
22 #include "AliRhoParameter.h"
23 #include "AliLog.h"
24 #include "AliEmcalParticle.h"
25 #include "AliMCEvent.h"
26 #include "AliGenPythiaEventHeader.h"
27 #include "AliAODMCHeader.h"
28 #include "AliMCEvent.h"
29 #include "AliAnalysisManager.h"
30 #include "AliJetContainer.h"
31
32 #include "AliAnalysisTaskEmcalDiJetBase.h"
33
34 ClassImp(AliAnalysisTaskEmcalDiJetBase)
35
36 //________________________________________________________________________
37 AliAnalysisTaskEmcalDiJetBase::AliAnalysisTaskEmcalDiJetBase() : 
38   AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalDiJetBase", kTRUE),
39   fDebug(kFALSE),
40   fJetCorrelationType(kCorrelateAll),
41   fJetFullChargedMatchingType(kFraction),
42   fTriggerClass(""),
43   fContainerCharged(1),
44   fContainerFull(0),
45   fContainerChargedMC(3),
46   fContainerFullMC(2),
47   fRhoType(0),
48   fRhoChVal(0),
49   fRhoFullVal(0),
50   fDoChargedCharged(kTRUE),
51   fDoFullCharged(kTRUE),
52   fDoFullFull(kFALSE),
53   fPtMinTriggerJet(10.),
54   fMinFractionShared(0.5),
55   fMatchingDone(kFALSE),
56   faFullFracIndex(0),
57   faFullFracIndexMC(0),
58   fhNEvents(0),
59   fHistTrialsSelEvents(0)
60 {
61   // Default constructor.
62
63   
64   SetMakeGeneralHistograms(kTRUE);
65
66 }
67
68 //________________________________________________________________________
69 AliAnalysisTaskEmcalDiJetBase::AliAnalysisTaskEmcalDiJetBase(const char *name) : 
70   AliAnalysisTaskEmcalJet(name, kTRUE),
71   fDebug(kFALSE),
72   fJetCorrelationType(kCorrelateAll),
73   fJetFullChargedMatchingType(kFraction),
74   fTriggerClass(""),
75   fContainerCharged(1),
76   fContainerFull(0),
77   fContainerChargedMC(3),
78   fContainerFullMC(2),
79   fRhoType(0),
80   fRhoChVal(0),
81   fRhoFullVal(0),
82   fDoChargedCharged(kTRUE),
83   fDoFullCharged(kTRUE),
84   fDoFullFull(kFALSE),
85   fPtMinTriggerJet(10.),
86   fMinFractionShared(0.5),
87   fMatchingDone(kFALSE),
88   faFullFracIndex(0),
89   faFullFracIndexMC(0),
90   fhNEvents(0),
91   fHistTrialsSelEvents(0)
92 {
93   // Standard constructor.
94
95   SetMakeGeneralHistograms(kTRUE);
96
97 }
98
99 //________________________________________________________________________
100 AliAnalysisTaskEmcalDiJetBase::~AliAnalysisTaskEmcalDiJetBase()
101 {
102   // Destructor.
103 }
104
105 //________________________________________________________________________
106 Bool_t AliAnalysisTaskEmcalDiJetBase::SelectEvent() {
107   //
108   // Decide if event should be selected for analysis
109   //
110
111   fhNEvents->Fill(3.5);
112
113   if(!fTriggerClass.IsNull()) {
114     //Check if requested trigger was fired
115     TString firedTrigClass = InputEvent()->GetFiredTriggerClasses();
116
117     if(!firedTrigClass.Contains(fTriggerClass))
118       return kFALSE;
119     else if(fTriggerClass.Contains("J1") && fTriggerClass.Contains("J2")) { //if events with J1&&J2 are requested
120       if(!firedTrigClass.Contains("J1") || !firedTrigClass.Contains("J2") ) //check if both are fired
121         return kFALSE;
122     }
123     else if(fTriggerClass.Contains("J1") && firedTrigClass.Contains("J2")) //if J2 is requested also add triggers which have J1&&J2. Reject if J1 is requested and J2 is fired
124       return kFALSE;
125   }
126   
127   fhNEvents->Fill(1.5);
128
129   fHistTrialsSelEvents->Fill(fPtHardBin, fNTrials);
130   
131   return kTRUE;
132
133 }
134
135 //________________________________________________________________________
136 void AliAnalysisTaskEmcalDiJetBase::UserCreateOutputObjects()
137 {
138   // Create user output.
139
140   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
141
142   Bool_t oldStatus = TH1::AddDirectoryStatus();
143   TH1::AddDirectory(kFALSE);
144
145   fhNEvents = new TH1F("fhNEvents","fhNEvents;selection;N_{evt}",5,0,5);
146   fOutput->Add(fhNEvents);
147
148   //Pythia info
149   fHistTrialsSelEvents = new TH1F("fHistTrialsSelEvents", "fHistTrialsSelEvents", 12, 0, 12);
150   fHistTrialsSelEvents->GetXaxis()->SetTitle("p_{T} hard bin");
151   fHistTrialsSelEvents->GetYaxis()->SetTitle("trials");
152   fOutput->Add(fHistTrialsSelEvents);
153
154   TH1::AddDirectory(oldStatus);
155
156   PostData(1, fOutput); // Post data for ALL output slots > 0 here.
157 }
158
159 //________________________________________________________________________
160 Bool_t AliAnalysisTaskEmcalDiJetBase::IsSameJet(const Int_t ijt, const Int_t ija, const Int_t type, const Bool_t isMC) {
161    //check if two jets are the same one
162
163    Bool_t bSame = kFALSE;
164
165    if(type<2 && ijt==ija)
166      bSame = kTRUE;
167    if(type==2) {
168
169      if(fJetFullChargedMatchingType==kFraction) {
170        if(isMC && faFullFracIndexMC.At(ijt)==ija)
171          bSame = kTRUE;
172        else if(!isMC && faFullFracIndex.At(ijt)==ija)
173          bSame = kTRUE;
174      }
175      else if(fJetFullChargedMatchingType==kGeo) {
176        Int_t contFull = fContainerFull;
177        Int_t contChar = fContainerCharged;
178        if(isMC) {
179          contFull = fContainerFullMC;
180          contChar = fContainerChargedMC;
181        }
182        AliEmcalJet *fullJet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ijt, contFull));
183        AliEmcalJet *chJet   = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ija, contChar));
184        AliEmcalJet *matchJet = fullJet->ClosestJet();
185        if(chJet==matchJet)
186          bSame = kTRUE;
187      }
188      else if(fJetFullChargedMatchingType==kNoMatching) {
189        return kFALSE;
190      }
191      else{
192        AliWarning(Form("%s: matching type unknown", GetName()));
193      }
194    }
195
196    return bSame;
197  }
198
199
200 //________________________________________________________________________
201 Double_t AliAnalysisTaskEmcalDiJetBase::GetJetPt(const AliEmcalJet *jet, const Int_t type) {
202
203   if(!jet) return -99;
204
205   if(type==0)
206     return jet->Pt() - fRhoFullVal * jet->Area();
207   else if(type==1)
208     return jet->Pt() - fRhoChVal * jet->Area();
209   else
210     return jet->Pt();
211   
212 }
213
214 //________________________________________________________________________
215 Double_t AliAnalysisTaskEmcalDiJetBase::GetZ(const AliVParticle *trk, const AliEmcalJet *jet)          const
216 {  
217   // Get Z of constituent trk
218
219   return GetZ(trk->Px(),trk->Py(),trk->Pz(),jet->Px(),jet->Py(),jet->Pz());
220 }
221
222 //________________________________________________________________________
223 Double_t AliAnalysisTaskEmcalDiJetBase::GetZ(const Double_t trkPx, const Double_t trkPy, const Double_t trkPz, const Double_t jetPx, const Double_t jetPy, const Double_t jetPz) const
224 {
225   // 
226   // Get the z of a constituent inside of a jet
227   //
228
229   Double_t pJetSq = jetPx*jetPx+jetPy*jetPy+jetPz*jetPz;
230
231   if(pJetSq>0.)
232     return (trkPx*jetPx+trkPy*jetPy+trkPz*jetPz)/pJetSq;
233   else {
234     AliWarning(Form("%s: strange, pjet*pjet seems to be zero pJetSq: %f",GetName(), pJetSq)); 
235     return 0;
236   }
237 }
238
239 //________________________________________________________________________
240 Double_t AliAnalysisTaskEmcalDiJetBase::GetDeltaR(const AliEmcalJet* jet1, const AliEmcalJet* jet2) const {
241   //
242   // Helper function to calculate the distance between two jets
243   //
244
245   Double_t dPhi = jet1->Phi() - jet2->Phi();
246   if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
247   if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
248   Double_t dEta = jet1->Eta() - jet2->Eta();
249   Double_t dR = TMath::Sqrt(dPhi*dPhi+dEta*dEta);
250   return dR;
251 }
252
253 //________________________________________________________________________
254 Double_t AliAnalysisTaskEmcalDiJetBase::GetDeltaPhi(const AliEmcalJet* jet1, const AliEmcalJet* jet2) {
255   //
256   // Calculate azimuthal angle between the axises of the jets
257   //
258
259   return GetDeltaPhi(jet1->Phi(),jet2->Phi());
260
261 }
262
263 //________________________________________________________________________
264 Double_t AliAnalysisTaskEmcalDiJetBase::GetDeltaPhi(Double_t phi1,Double_t phi2) {
265   //
266   // Calculate azimuthal angle between the axises of the jets
267   //
268
269   Double_t dPhi = phi1-phi2;
270
271   if(dPhi <-0.5*TMath::Pi())  dPhi += TMath::TwoPi();
272   if(dPhi > 1.5*TMath::Pi())  dPhi -= TMath::TwoPi();
273
274   return dPhi;
275 }
276
277 //________________________________________________________________________
278 Double_t AliAnalysisTaskEmcalDiJetBase::GetFractionSharedPt(const AliEmcalJet *jetFull, const AliEmcalJet *jetCharged) const
279 {
280   //
281   // Get fraction of shared pT between matched full and charged jet
282   // Uses charged jet pT as baseline: fraction = \Sum_{const,full jet} pT,const,i / pT,jet,ch
283   //
284
285   Double_t fraction = 0.;
286   Double_t jetPtCh = jetCharged->Pt();
287  
288   if(jetPtCh>0) {
289
290     AliDebug(11,Form("%s: nConstituents: %d, ch: %d  chne: %d ne: %d",GetName(),jetFull->GetNumberOfConstituents(),jetCharged->GetNumberOfTracks(),jetFull->GetNumberOfTracks(),jetFull->GetNumberOfClusters()));
291     
292     Double_t sumPt = 0.;
293     AliVParticle *vpf = 0x0;
294     Int_t iFound = 0;
295     for(Int_t icc=0; icc<jetCharged->GetNumberOfTracks(); icc++) {
296       Int_t idx = (Int_t)jetCharged->TrackAt(icc);
297       iFound = 0;
298       for(Int_t icf=0; icf<jetFull->GetNumberOfTracks(); icf++) {
299         if(idx == jetFull->TrackAt(icf) && iFound==0 ) {
300           iFound=1;
301           vpf = static_cast<AliVParticle*>(jetFull->TrackAt(icf, fTracks));
302           if(!vpf) continue;
303           if(vpf->Charge()!=0)
304             sumPt += vpf->Pt();
305           continue;
306         }
307       }
308     }
309     
310     fraction = sumPt/jetPtCh;
311   }
312
313   AliDebug(11,Form("%s: charged shared fraction: %.2f",GetName(),fraction));
314
315   return fraction;
316
317 }
318
319 //________________________________________________________________________
320 void AliAnalysisTaskEmcalDiJetBase::MatchJetsGeo(Int_t cFull, Int_t cCharged,
321                                                  Int_t iDebug, Float_t maxDist, Int_t type) {
322
323   //
324   // Match the full jets to the corresponding charged jets
325   // Translation of AliAnalysisHelperJetTasks::GetClosestJets to AliEmcalJet objects
326   // type: 0 = use acceptance cuts of container  1 = allow 0.1 one more for cCharged(MC) in eta 2 = allow 0.1 more in eta and phi for cCharged(MC)
327
328   const int kMode = 3;
329
330   const Int_t nChJets   = GetNJets(cCharged);
331   const Int_t nFullJets = GetNJets(cFull);
332
333   if(nFullJets==0 || nChJets==0) return;
334
335   TArrayI faChNeMatchIndex;
336   faChNeMatchIndex.Set(nChJets+1);
337   faChNeMatchIndex.Reset(-1);
338
339   TArrayI faChMatchIndex;
340   faChMatchIndex.Set(nFullJets+1);
341   faChMatchIndex.Reset(-1);
342
343   static TArrayS iFlag(nChJets*nFullJets);
344   if(iFlag.GetSize()<(nChJets*nFullJets)){
345     iFlag.Set(nChJets*nFullJets+1);
346   }
347   iFlag.Reset(0);
348
349   AliJetContainer *contCh = GetJetContainer(cCharged);
350
351   // find the closest distance to the full jet
352   for(int ifu = 0;ifu<nFullJets;ifu++){
353
354     AliEmcalJet *fullJet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ifu, cFull));
355     if(!fullJet) continue;
356
357     Float_t dist = maxDist;
358     
359     for(int ich = 0;ich<nChJets;ich++){
360       AliEmcalJet *chJet = 0x0;
361       if(type==0)
362         chJet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ich, cCharged));
363       else {
364         chJet = static_cast<AliEmcalJet*>(GetJetFromArray(ich, cCharged));
365         if(!chJet) continue;
366         if(type>0) {
367           if(chJet->Eta()<(contCh->GetJetEtaMin()-0.1) || chJet->Eta()>(contCh->GetJetEtaMax()+0.1))
368             continue;
369           if(type==2) {
370             if(chJet->Phi()<(contCh->GetJetPhiMin()-0.1) || chJet->Phi()>(contCh->GetJetPhiMax()+0.1))
371               continue;
372           }
373         }
374       }
375       if(!chJet)
376         continue;
377
378       Double_t frac = GetFractionSharedPt(fullJet,chJet);
379       Double_t dR = GetDeltaR(fullJet,chJet);
380       if(dR<dist){
381         faChMatchIndex[ifu]=ich;
382         dist = dR;
383       }
384       if(iDebug>10) Printf("Distance (%d)--(%d) %3.3f frac:%.2f",ifu,ich,dR,frac);
385     }
386     if(faChMatchIndex[ifu]>=0) iFlag[ifu*nChJets+faChMatchIndex[ifu]]+=1;//ich closest to ifu
387     if(iDebug>10) Printf("Full Distance (%d)--(%d) %3.3f flag[%d] = %d",ifu,faChMatchIndex[ifu],dist,ifu*nChJets+faChMatchIndex[ifu],iFlag[ifu*nChJets+faChMatchIndex[ifu]]);
388     
389     // reset...
390     faChMatchIndex[ifu]=-1;
391
392            
393   }//full jet loop
394
395
396   // other way around
397   for(int ich = 0;ich<nChJets;ich++){
398     AliEmcalJet *chJet = 0x0;
399     if(type==0)
400       chJet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ich, cCharged));
401     else {
402       chJet = static_cast<AliEmcalJet*>(GetJetFromArray(ich, cCharged));
403       if(!chJet) continue;;
404       if(type>0) {
405         if(chJet->Eta()<(contCh->GetJetEtaMin()-0.1) || chJet->Eta()>(contCh->GetJetEtaMax()+0.1))
406           continue;
407         if(type==2) {
408           if(chJet->Phi()<(contCh->GetJetPhiMin()-0.1) || chJet->Phi()>(contCh->GetJetPhiMax()+0.1))
409             continue;
410         }
411       }
412     }
413     if(!chJet)
414       continue;
415
416     Float_t dist = maxDist;
417     for(int ifu = 0;ifu<nFullJets;ifu++){
418       AliEmcalJet *fullJet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ifu, cFull));
419       if(!fullJet)
420         continue;
421
422       Double_t dR = GetDeltaR(fullJet,chJet); 
423       if(dR<dist){
424         faChNeMatchIndex[ich]=ifu;
425         dist = dR;
426       }   
427     }
428     if(faChNeMatchIndex[ich]>=0) iFlag[faChNeMatchIndex[ich]*nChJets+ich]+=2;//ifu closest to ich
429     if(iDebug>10) Printf("Other way Distance (%d)--(%d) %3.3f flag[%d] = %d",faChNeMatchIndex[ich],ich,dist,faChNeMatchIndex[ich]*nChJets+ich,iFlag[faChNeMatchIndex[ich]*nChJets+ich]);
430
431     // reset...
432     faChNeMatchIndex[ich]=-1;
433     
434   }
435   
436   
437   // check for "true" correlations
438   for(int ifu = 0;ifu<nFullJets;ifu++){
439     for(int ich = 0;ich<nChJets;ich++){
440       AliDebug(11,Form("%s: Flag[%d][%d] %d ",GetName(),ifu,ich,iFlag[ifu*nChJets+ich]));
441       
442       if(kMode==3){
443         // we have a uniqe correlation
444         if(iFlag[ifu*nChJets+ich]==3){
445
446           AliEmcalJet *chJet = static_cast<AliEmcalJet*>(GetJetFromArray(ich, cCharged));
447           AliEmcalJet *fullJet = static_cast<AliEmcalJet*>(GetJetFromArray(ifu, cFull));
448           Double_t dR = GetDeltaR(fullJet,chJet); 
449
450           AliDebug(11,Form("closest jets %d  %d  dR =  %f",ich,ifu,dR));
451
452           chJet->SetClosestJet(fullJet,dR);
453           fullJet->SetClosestJet(chJet,dR);
454
455         }
456       }
457     }
458   }
459
460   fMatchingDone = kTRUE;
461   
462 }
463
464 //________________________________________________________________________
465 void AliAnalysisTaskEmcalDiJetBase::SetChargedFractionIndex() {
466
467   // take each full jet and set the index of the charged jet with the largest shared charged fraction
468
469   const Int_t nJetsCh   = GetNJets(fContainerCharged);
470   const Int_t nJetsFull = GetNJets(fContainerFull);
471   faFullFracIndex.Set(nJetsFull+1);
472   faFullFracIndex.Reset(-1);
473
474   AliJetContainer  *cont = GetJetContainer(fContainerFull); 
475   Float_t radius =  cont->GetJetRadius();
476
477   for(Int_t ifu = 0; ifu<nJetsFull; ifu++) {
478     Double_t frac = 0.;
479     Double_t dist = 10.;
480     AliEmcalJet *jetFull = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ifu, fContainerFull));
481     if(!jetFull) {
482       faFullFracIndex.SetAt(-1,ifu);
483       continue;
484     }
485
486     for(Int_t ich = 0; ich<nJetsCh; ich++) {
487       AliEmcalJet *jetCh = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ich, fContainerCharged));
488       if(!jetCh)
489         continue;
490       Double_t tmpFrac = GetFractionSharedPt(jetFull,jetCh);
491       dist = GetDeltaR(jetFull,jetCh);
492       if(tmpFrac>frac && dist<radius) {
493         frac = tmpFrac;
494         faFullFracIndex.SetAt(ich,ifu);
495       }
496     }
497
498   }
499
500 }
501
502 //________________________________________________________________________
503 void AliAnalysisTaskEmcalDiJetBase::SetChargedFractionIndexMC() {
504
505   // take each full jet and set the index of the charged jet with the largest shared charged fraction
506
507   const Int_t nJetsCh   = GetNJets(fContainerChargedMC);
508   const Int_t nJetsFull = GetNJets(fContainerFullMC);
509   faFullFracIndexMC.Set(nJetsFull);
510   faFullFracIndexMC.Reset(-1);
511
512   AliJetContainer  *cont = GetJetContainer(fContainerFullMC); 
513   Float_t radius =  cont->GetJetRadius();
514
515   for(Int_t ifu = 0; ifu<nJetsFull; ifu++) {
516     Double_t frac = 0.;
517     Double_t dist = 10.;
518     AliEmcalJet *jetFull = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ifu, fContainerFullMC));
519     if(!jetFull) {
520       faFullFracIndexMC.SetAt(-1,ifu);
521       continue;
522     }
523
524     for(Int_t ich = 0; ich<nJetsCh; ich++) {
525       AliEmcalJet *jetCh = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ich, fContainerChargedMC));
526       if(!jetCh)
527         continue;
528       Double_t tmpFrac = GetFractionSharedPt(jetFull,jetCh);
529       dist = GetDeltaR(jetFull,jetCh);
530       if(tmpFrac>frac && dist<radius) {
531         frac = tmpFrac;
532         faFullFracIndexMC.SetAt(ich,ifu);
533       }
534     }
535
536   }
537
538 }
539
540 //_______________________________________________________________________
541 AliEmcalJet* AliAnalysisTaskEmcalDiJetBase::GetLeadingJetOppositeHemisphere(const Int_t type, const Int_t typea, const AliEmcalJet *jetTrig) {
542
543   // Get leading jet in opposite hemisphere from trigger jet
544   // type = correlation type
545   // typea = container of associated jets
546
547   Int_t nJetsAssoc = GetNJets(typea);
548   Double_t ptLead = -999;
549   Int_t    iJetLead = -1;
550   for(Int_t ija=0; ija<nJetsAssoc; ija++) {
551
552     AliEmcalJet *jetAssoc = NULL;
553     if(type==0) {
554       jetAssoc = static_cast<AliEmcalJet*>(GetJetFromArray(ija, typea));
555       if(TMath::Abs(jetAssoc->Eta())>0.5)
556         jetAssoc = NULL;
557     }
558     else
559       jetAssoc = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ija, typea));
560
561     if(!jetAssoc)
562       continue;
563
564     Double_t dPhi = GetDeltaPhi(jetTrig,jetAssoc);
565     Double_t phiMin = 0.5*TMath::Pi();
566     Double_t phiMax = 1.5*TMath::Pi();
567     if(dPhi<phiMin || dPhi>phiMax)
568       continue;
569
570     Double_t jetAssocPt = GetJetPt(jetAssoc,typea);
571
572     if(jetAssocPt>ptLead) {
573       ptLead = jetAssocPt;
574       iJetLead = ija;
575     }
576
577   }
578
579   AliEmcalJet *jetAssocLead = NULL;
580   if(iJetLead>-1)
581     jetAssocLead = static_cast<AliEmcalJet*>(GetJetFromArray(iJetLead, typea));
582
583   return jetAssocLead;
584
585 }
586
587 //_______________________________________________________________________
588 AliEmcalJet* AliAnalysisTaskEmcalDiJetBase::GetSecondLeadingJetOppositeHemisphere(const Int_t type, const Int_t typea, const AliEmcalJet *jetTrig) {
589
590   // Get leading jet in opposite hemisphere from trigger jet
591   // type = correlation type
592   // typea = container of associated jets
593
594   Int_t nJetsAssoc = GetNJets(typea);
595   Double_t ptLead = -999;
596   Int_t    iJetLead = -1;
597   Int_t    iJetLead2 = -1;
598   for(Int_t ija=0; ija<nJetsAssoc; ija++) {
599
600     AliEmcalJet *jetAssoc = NULL;
601     if(type==0) {
602       jetAssoc = static_cast<AliEmcalJet*>(GetJetFromArray(ija, typea));
603       if(TMath::Abs(jetAssoc->Eta())>0.5)
604         jetAssoc = NULL;
605     }
606     else
607       jetAssoc = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ija, typea));
608
609     if(!jetAssoc)
610       continue;
611
612     Double_t dPhi = GetDeltaPhi(jetTrig,jetAssoc);
613     Double_t phiMin = 0.5*TMath::Pi();
614     Double_t phiMax = 1.5*TMath::Pi();
615     if(dPhi<phiMin || dPhi>phiMax)
616       continue;
617
618     Double_t jetAssocPt = GetJetPt(jetAssoc,typea);
619
620     if(jetAssocPt>ptLead) {
621       iJetLead2 = iJetLead;
622       ptLead = jetAssocPt;
623       iJetLead = ija;
624     }
625
626   }
627
628   AliEmcalJet *jetAssocLead2 = NULL;
629   if(iJetLead2>-1)
630     jetAssocLead2 = static_cast<AliEmcalJet*>(GetJetFromArray(iJetLead2, typea));
631
632   return jetAssocLead2;
633
634 }
635
636 //________________________________________________________________________
637 Bool_t AliAnalysisTaskEmcalDiJetBase::RetrieveEventObjects() {
638   //
639   // Retrieve objects from event.
640   //
641
642   if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
643     return kFALSE;
644
645   if(fRhoType==0) {
646     fRhoFullVal = 0.;
647     fRhoChVal = 0.;
648   }
649   if(fRhoType==1) {
650     fRhoFullVal = GetRhoVal(fContainerFull);
651     fRhoChVal = GetRhoVal(fContainerCharged);
652   }
653
654   return kTRUE;
655 }
656
657 //_______________________________________________________________________
658 void AliAnalysisTaskEmcalDiJetBase::Terminate(Option_t *) 
659 {
660   // Called once at the end of the analysis.
661 }