PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetTagger.cxx

   1 //
   2 // Jet tagger 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 #include "AliParticleContainer.h"
  32
  33 #include "AliAODEvent.h"
  34 #include "AliESDEvent.h"
  35
  36 #include "AliAnalysisTaskEmcalJetTagger.h"
  37
  38 ClassImp(AliAnalysisTaskEmcalJetTagger)
  39
  40 //________________________________________________________________________
  41 AliAnalysisTaskEmcalJetTagger::AliAnalysisTaskEmcalJetTagger() :
  42   AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetTagger", kTRUE),
  43   fJetTaggingType(kTag),
  44   fJetTaggingMethod(kGeo),
  45   fContainerBase(0),
  46   fContainerTag(1),
  47   fMinFractionShared(0),
  48   fMatchingDone(0),
  49   fh3PtJet1VsDeltaEtaDeltaPhi(0),
  50   fh2PtJet1VsDeltaR(0),
  51   fh2PtJet2VsFraction(0),
  52   fh2PtJet1VsLeadPtAllSel(0),
  53   fh2PtJet1VsLeadPtTagged(0),
  54   fh2PtJet1VsPtJet2(0),
  55   fh2PtJet2VsRelPt(0),
  56   fh3PtJetDEtaDPhiConst(0),
  57   fh2PtJetDRConst(0),
  58   fh3PtJetAreaDRConst(0)
  59 {
  60   // Default constructor.
  61
  62   fh3PtJet1VsDeltaEtaDeltaPhi  = new TH3F*[fNcentBins];
  63   fh2PtJet1VsDeltaR            = new TH2F*[fNcentBins];
  64   fh2PtJet2VsFraction          = new TH2F*[fNcentBins];
  65   fh2PtJet1VsLeadPtAllSel      = new TH2F*[fNcentBins];
  66   fh2PtJet1VsLeadPtTagged      = new TH2F*[fNcentBins];
  67   fh2PtJet1VsPtJet2            = new TH2F*[fNcentBins];
  68   fh2PtJet2VsRelPt             = new TH2F*[fNcentBins];
  69
  70   for (Int_t i = 0; i < fNcentBins; i++) {
  71     fh3PtJet1VsDeltaEtaDeltaPhi[i] = 0;
  72     fh2PtJet1VsDeltaR[i]           = 0;
  73     fh2PtJet2VsFraction[i]         = 0;
  74     fh2PtJet1VsLeadPtAllSel[i]     = 0;
  75     fh2PtJet1VsLeadPtTagged[i]     = 0;
  76     fh2PtJet1VsPtJet2[i]           = 0;
  77     fh2PtJet2VsRelPt[i]            = 0;
  78   }
  79
  80   SetMakeGeneralHistograms(kTRUE);
  81 }
  82
  83 //________________________________________________________________________
  84 AliAnalysisTaskEmcalJetTagger::AliAnalysisTaskEmcalJetTagger(const char *name) :
  85   AliAnalysisTaskEmcalJet(name, kTRUE),
  86   fJetTaggingType(kTag),
  87   fJetTaggingMethod(kGeo),
  88   fContainerBase(0),
  89   fContainerTag(1),
  90   fMinFractionShared(0),
  91   fMatchingDone(0),
  92   fh3PtJet1VsDeltaEtaDeltaPhi(0),
  93   fh2PtJet1VsDeltaR(0),
  94   fh2PtJet2VsFraction(0),
  95   fh2PtJet1VsLeadPtAllSel(0),
  96   fh2PtJet1VsLeadPtTagged(0),
  97   fh2PtJet1VsPtJet2(0),
  98   fh2PtJet2VsRelPt(0),
  99   fh3PtJetDEtaDPhiConst(0),
 100   fh2PtJetDRConst(0),
 101   fh3PtJetAreaDRConst(0)
 102 {
 103   // Standard constructor.
 104
 105   fh3PtJet1VsDeltaEtaDeltaPhi = new TH3F*[fNcentBins];
 106   fh2PtJet1VsDeltaR           = new TH2F*[fNcentBins];
 107   fh2PtJet2VsFraction         = new TH2F*[fNcentBins];
 108   fh2PtJet1VsLeadPtAllSel     = new TH2F*[fNcentBins];
 109   fh2PtJet1VsLeadPtTagged     = new TH2F*[fNcentBins];
 110   fh2PtJet1VsPtJet2           = new TH2F*[fNcentBins];
 111   fh2PtJet2VsRelPt            = new TH2F*[fNcentBins];
 112
 113   for (Int_t i = 0; i < fNcentBins; i++) {
 114     fh3PtJet1VsDeltaEtaDeltaPhi[i] = 0;
 115     fh2PtJet1VsDeltaR[i]           = 0;
 116     fh2PtJet2VsFraction[i]         = 0;
 117     fh2PtJet1VsLeadPtAllSel[i]     = 0;
 118     fh2PtJet1VsLeadPtTagged[i]     = 0;
 119     fh2PtJet1VsPtJet2[i]           = 0;
 120     fh2PtJet2VsRelPt[i]            = 0;
 121   }
 122
 123   SetMakeGeneralHistograms(kTRUE);
 124 }
 125
 126 //________________________________________________________________________
 127 AliAnalysisTaskEmcalJetTagger::~AliAnalysisTaskEmcalJetTagger()
 128 {
 129   // Destructor.
 130 }
 131
 132 //________________________________________________________________________
 133 void AliAnalysisTaskEmcalJetTagger::UserCreateOutputObjects()
 134 {
 135   // Create user output.
 136
 137   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 138
 139   Bool_t oldStatus = TH1::AddDirectoryStatus();
 140   TH1::AddDirectory(kFALSE);
 141
 142   const Int_t nBinsPt          = 250;
 143   const Int_t nBinsDPhi        = 100;
 144   const Int_t nBinsDEta        = 100;
 145   const Int_t nBinsDR          = 50;
 146   const Int_t nBinsFraction    = 101;
 147
 148   const Double_t minPt       = -50.;
 149   const Double_t maxPt       = 200.;
 150   const Double_t minDPhi     = -0.5;
 151   const Double_t maxDPhi     =  0.5;
 152   const Double_t minDEta     = -0.5;
 153   const Double_t maxDEta     =  0.5;
 154   const Double_t minDR       =  0.;
 155   const Double_t maxDR       =  0.5;
 156   const Double_t minFraction =  -0.005;
 157   const Double_t maxFraction =  1.005;
 158
 159   TString histName = "";
 160   TString histTitle = "";
 161   for (Int_t i = 0; i < fNcentBins; i++) {
 162
 163     histName = TString::Format("fh3PtJet1VsDeltaEtaDeltaPhi_%d",i);
 164     histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{#Delta#eta};#it{#Delta#varphi}",histName.Data());
 165     fh3PtJet1VsDeltaEtaDeltaPhi[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsDEta,minDEta,maxDEta,nBinsDPhi,minDPhi,maxDPhi);
 166     fOutput->Add(fh3PtJet1VsDeltaEtaDeltaPhi[i]);
 167
 168     histName = TString::Format("fh2PtJet1VsDeltaR_%d",i);
 169     histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{#Delta R}",histName.Data());
 170     fh2PtJet1VsDeltaR[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsDR,minDR,maxDR);
 171     fOutput->Add(fh2PtJet1VsDeltaR[i]);
 172
 173     histName = TString::Format("fh2PtJet2VsFraction_%d",i);
 174     histTitle = TString::Format("%s;#it{p}_{T,jet2};#it{f}_{shared}",histName.Data());
 175     fh2PtJet2VsFraction[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsFraction,minFraction,maxFraction);
 176     fOutput->Add(fh2PtJet2VsFraction[i]);
 177
 178     histName = TString::Format("fh2PtJet1VsLeadPtAllSel_%d",i);
 179     histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{p}_{T,lead trk}",histName.Data());
 180     fh2PtJet1VsLeadPtAllSel[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,20,0.,20.);
 181     fOutput->Add(fh2PtJet1VsLeadPtAllSel[i]);
 182
 183     histName = TString::Format("fh2PtJet1VsLeadPtTagged_%d",i);
 184     histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{p}_{T,lead trk}",histName.Data());
 185     fh2PtJet1VsLeadPtTagged[i] =  new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,20,0.,20.);
 186     fOutput->Add(fh2PtJet1VsLeadPtTagged[i]);
 187
 188     histName = TString::Format("fh2PtJet1VsPtJet2_%d",i);
 189     histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{p}_{T,jet2}",histName.Data());
 190     fh2PtJet1VsPtJet2[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt);
 191     fOutput->Add(fh2PtJet1VsPtJet2[i]);
 192
 193     histName = TString::Format("fh2PtJet2VsRelPt_%d",i);
 194     histTitle = TString::Format("%s;#it{p}_{T,jet2};(#it{p}_{T,jet2}-#it{p}_{T,jet1})/#it{p}_{T,jet1}",histName.Data());
 195     fh2PtJet2VsRelPt[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,241,-2.41,2.41);
 196     fOutput->Add(fh2PtJet2VsRelPt[i]);
 197   }
 198
 199   fh3PtJetDEtaDPhiConst = new TH3F("fh3PtJetDEtaDPhiConst","fh3PtJetDEtaDPhiConst;pT;#Delta #eta;#Delta #varphi",nBinsPt,minPt,maxPt,nBinsDEta,-1.,1.,nBinsDPhi,-1.,1.);
 200   fOutput->Add(fh3PtJetDEtaDPhiConst);
 201
 202   fh2PtJetDRConst = new TH2F("fh2PtJetDRConst","fh2PtJetDRConst;pT;#Delta R",nBinsPt,minPt,maxPt,100,0.,1.);
 203   fOutput->Add(fh2PtJetDRConst);
 204
 205   fh3PtJetAreaDRConst = new TH3F("fh3PtJetAreaDRConst","fh3PtJetAreaDRConst;pT;A;#Delta R",nBinsPt,minPt,maxPt,100,0.,1.,100,0.,1.);
 206   fOutput->Add(fh3PtJetAreaDRConst);
 207
 208   // =========== Switch on Sumw2 for all histos ===========
 209   for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
 210     TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
 211     if (h1){
 212       h1->Sumw2();
 213       continue;
 214     }
 215     THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
 216     if(hn)hn->Sumw2();
 217   }
 218
 219   TH1::AddDirectory(oldStatus);
 220
 221   PostData(1, fOutput); // Post data for ALL output slots > 0 here.
 222 }
 223
 224 //________________________________________________________________________
 225 Bool_t AliAnalysisTaskEmcalJetTagger::Run()
 226 {
 227   // Run analysis code here, if needed. It will be executed before FillHistograms().
 228
 229   MatchJetsGeo(fContainerBase,fContainerTag,0,0.3,3);
 230
 231   //  if(fJetTaggingMethod==kFraction)
 232
 233   return kTRUE;
 234 }
 235
 236 //________________________________________________________________________
 237 Bool_t AliAnalysisTaskEmcalJetTagger::FillHistograms()
 238 {
 239   // Fill histograms.
 240
 241   AliEmcalJet *jet1 = NULL;
 242   AliJetContainer *jetCont = GetJetContainer(fContainerBase);
 243   if(!jetCont) return kFALSE;
 244   jetCont->ResetCurrentID();
 245   while((jet1 = jetCont->GetNextAcceptJet())) {
 246     Double_t ptJet1 =  jet1->Pt() - jetCont->GetRhoVal()*jet1->Area();
 247     fh2PtJet1VsLeadPtAllSel[fCentBin]->Fill(ptJet1,jet1->MaxTrackPt());
 248
 249     //fill histo with angle between jet axis and constituents
 250     for(Int_t icc=0; icc<jet1->GetNumberOfTracks(); icc++) {
 251       AliVParticle *vp = static_cast<AliVParticle*>(jet1->TrackAt(icc, jetCont->GetParticleContainer()->GetArray()));//fTracks));
 252       if(!vp) continue;
 253       Double_t dEta = jet1->Eta()-vp->Eta();
 254       Double_t dPhi = jet1->Phi()-vp->Phi();
 255       if(dPhi<TMath::Pi()) dPhi+=TMath::TwoPi();
 256       if(dPhi>TMath::Pi()) dPhi-=TMath::TwoPi();
 257       fh3PtJetDEtaDPhiConst->Fill(ptJet1,dEta,dPhi);
 258
 259       Double_t dR = TMath::Sqrt(dPhi*dPhi+dEta*dEta);
 260       fh2PtJetDRConst->Fill(ptJet1,dR);
 261       fh3PtJetAreaDRConst->Fill(ptJet1,jet1->Area(),dR);
 262     }
 263
 264     if(jet1->GetTagStatus()<1 && fJetTaggingType==kTag)
 265       continue;
 266
 267     AliEmcalJet *jet2 = NULL;
 268     if(fJetTaggingType==kTag)     jet2 = jet1->GetTaggedJet();
 269     if(fJetTaggingType==kClosest) jet2 = jet1->ClosestJet();
 270     if(!jet2) continue;
 271
 272     Double_t ptJet2 =  jet2->Pt() - GetRhoVal(fContainerTag)*jet2->Area();
 273     Double_t fraction = jetCont->GetFractionSharedPt(jet1);
 274     fh2PtJet2VsFraction[fCentBin]->Fill(ptJet2,fraction);
 275
 276     if(fraction<fMinFractionShared && fJetTaggingType==kClosest)
 277       continue;
 278
 279     fh2PtJet1VsLeadPtTagged[fCentBin]->Fill(ptJet1,jet1->MaxTrackPt());
 280     fh2PtJet1VsPtJet2[fCentBin]->Fill(ptJet1,ptJet2);
 281     if(ptJet2>0.) fh2PtJet2VsRelPt[fCentBin]->Fill(ptJet2,(ptJet1-ptJet2)/ptJet2);
 282
 283     Double_t dPhi = GetDeltaPhi(jet1->Phi(),jet2->Phi());
 284     if(dPhi>TMath::Pi())
 285       dPhi -= TMath::TwoPi();
 286     if(dPhi<(-1.*TMath::Pi()))
 287       dPhi += TMath::TwoPi();
 288
 289     fh3PtJet1VsDeltaEtaDeltaPhi[fCentBin]->Fill(ptJet1,jet1->Eta()-jet2->Eta(),dPhi);
 290     fh2PtJet1VsDeltaR[fCentBin]->Fill(ptJet1,GetDeltaR(jet1,jet2));
 291   }
 292   return kTRUE;
 293 }
 294
 295 //________________________________________________________________________
 296 void AliAnalysisTaskEmcalJetTagger::ResetTagging(const Int_t c) {
 297
 298   //Reset tagging of container c
 299
 300   for(int i = 0;i<GetNJets(c);i++){
 301     AliEmcalJet *jet = static_cast<AliEmcalJet*>(GetJetFromArray(i, c));
 302     if(fJetTaggingType==kClosest)
 303       jet->ResetMatching();
 304     else if(fJetTaggingType==kTag) {
 305       jet->SetTaggedJet(0x0);
 306       jet->SetTagStatus(-1);
 307     }
 308   }
 309 }
 310
 311 //________________________________________________________________________
 312 void AliAnalysisTaskEmcalJetTagger::MatchJetsGeo(Int_t c1, Int_t c2,
 313                                                  Int_t iDebug, Float_t maxDist, Int_t type, Bool_t bReset) {
 314
 315   //
 316   // Match the full jets to the corresponding charged jets
 317   // Translation of AliAnalysisHelperJetTasks::GetClosestJets to AliEmcalJet objects
 318   // type:
 319   //         0 = use acceptance cuts of container
 320   //         1 = allow 0.1 one more for c2 in eta
 321   //         2 = allow 0.1 more in eta and phi for c2
 322   //         3 = allow 0.1 in eta and phi for both containers
 323
 324   if(c1<0) c1 = fContainerBase;
 325   if(c2<0) c2 = fContainerTag;
 326
 327   const Int_t nJets1 = GetNJets(c1);
 328   const Int_t nJets2 = GetNJets(c2);
 329
 330   if(nJets1==0 || nJets2==0) return;
 331
 332   if(bReset) {
 333     ResetTagging(c1);
 334     ResetTagging(c2);
 335   }
 336   fMatchingDone = kFALSE;
 337
 338   TArrayI faMatchIndex1;
 339   faMatchIndex1.Set(nJets2+1);
 340   faMatchIndex1.Reset(-1);
 341
 342   TArrayI faMatchIndex2;
 343   faMatchIndex2.Set(nJets1+1);
 344   faMatchIndex2.Reset(-1);
 345
 346   static TArrayS iFlag(nJets1*nJets2);
 347   if(iFlag.GetSize()<(nJets1*nJets2)){
 348     iFlag.Set(nJets1*nJets2+1);
 349   }
 350   iFlag.Reset(0);
 351
 352
 353   AliJetContainer *cont1 = GetJetContainer(c1);
 354   AliJetContainer *cont2 = GetJetContainer(c2);
 355   if(type==1)
 356     cont2->SetJetEtaLimits(cont2->GetJetEtaMin()-0.1,cont2->GetJetEtaMax()+0.1);
 357   else if(type==2) {
 358     cont2->SetJetEtaLimits(cont2->GetJetEtaMin()-0.1,cont2->GetJetEtaMax()+0.1);
 359     cont2->SetJetPhiLimits(cont2->GetJetPhiMin()-0.1,cont2->GetJetPhiMax()+0.1);
 360   }
 361   else if(type==3) {
 362     cont1->SetJetEtaLimits(cont1->GetJetEtaMin()-0.1,cont1->GetJetEtaMax()+0.1);
 363     cont1->SetJetPhiLimits(cont1->GetJetPhiMin()-0.1,cont1->GetJetPhiMax()+0.1);
 364     cont2->SetJetEtaLimits(cont2->GetJetEtaMin()-0.1,cont2->GetJetEtaMax()+0.1);
 365     cont2->SetJetPhiLimits(cont2->GetJetPhiMin()-0.1,cont2->GetJetPhiMax()+0.1);
 366   }
 367
 368   // find the closest distance to the full jet
 369   for(int i = 0;i<nJets1;i++){
 370
 371     AliEmcalJet *jet1 = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(i, c1));
 372     if(!jet1) continue;
 373
 374     Float_t dist = maxDist;
 375
 376     for(int j = 0;j <nJets2; j++){
 377       AliEmcalJet *jet2 = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(j, c2));
 378       if(!jet2) continue;
 379
 380       Double_t dR = jet1->DeltaR(jet2);//GetDeltaR(jet1,jet2);
 381       if(dR<dist && dR<maxDist){
 382         faMatchIndex2[i]=j;
 383         dist = dR;
 384       }
 385     }//j jet loop
 386     if(faMatchIndex2[i]>=0) {
 387       iFlag[i*nJets2+faMatchIndex2[i]]+=1;//j closest to i
 388       if(iDebug>10) Printf("Full Distance (%d)--(%d) %3.3f flag[%d] = %d",i,faMatchIndex2[i],dist,i*nJets2+faMatchIndex2[i],iFlag[i*nJets2+faMatchIndex2[i]]);
 389     }
 390
 391   }//i jet loop
 392
 393
 394   // other way around
 395   for(int j = 0;j<nJets2;j++){
 396     AliEmcalJet *jet2 = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(j, c2));
 397     if(!jet2)
 398       continue;
 399
 400     Float_t dist = maxDist;
 401     for(int i = 0;i<nJets1;i++){
 402       AliEmcalJet *jet1 = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(i, c1));
 403       if(!jet1) continue;
 404
 405       Double_t dR = jet1->DeltaR(jet2);//GetDeltaR(jet1,jet2);
 406       if(dR<dist && dR<maxDist){
 407         faMatchIndex1[j]=i;
 408         dist = dR;
 409       }
 410     }
 411     if(faMatchIndex1[j]>=0) {
 412       iFlag[faMatchIndex1[j]*nJets2+j]+=2;//i closest to j
 413       if(iDebug>10) Printf("Other way Distance (%d)--(%d) %3.3f flag[%d] = %d",faMatchIndex1[j],j,dist,faMatchIndex1[j]*nJets2+j,iFlag[faMatchIndex1[j]*nJets2+j]);
 414     }
 415
 416   }
 417
 418   // check for "true" correlations
 419   for(int i = 0;i<nJets1;i++){
 420     AliEmcalJet *jet1 = static_cast<AliEmcalJet*>(GetJetFromArray(i, c1));
 421     for(int j = 0;j<nJets2;j++){
 422       AliEmcalJet *jet2 = static_cast<AliEmcalJet*>(GetJetFromArray(j, c2));
 423       AliDebug(11,Form("%s: Flag[%d][%d] %d ",GetName(),i,j,iFlag[i*nJets2+j]));
 424
 425       // we have a uniqe correlation
 426       if(iFlag[i*nJets2+j]==3){
 427
 428         Double_t dR = jet1->DeltaR(jet2);//GetDeltaR(jet1,jet2);
 429         if(iDebug>1) Printf("closest jets %d  %d  dR =  %f",j,i,dR);
 430
 431         if(fJetTaggingType==kTag) {
 432           jet1->SetTaggedJet(jet2);
 433           jet1->SetTagStatus(1);
 434
 435           jet2->SetTaggedJet(jet1);
 436           jet2->SetTagStatus(1);
 437         }
 438         else if(fJetTaggingType==kClosest) {
 439           jet1->SetClosestJet(jet2,dR);
 440           jet2->SetClosestJet(jet1,dR);
 441         }
 442       }
 443     }
 444   }
 445   fMatchingDone = kTRUE;
 446 }
 447
 448 //________________________________________________________________________
 449 Double_t AliAnalysisTaskEmcalJetTagger::GetDeltaR(const AliEmcalJet* jet1, const AliEmcalJet* jet2) const {
 450   //
 451   // Helper function to calculate the distance between two jets
 452   //
 453
 454   Double_t dPhi = jet1->Phi() - jet2->Phi();
 455   if(dPhi>TMath::Pi())
 456     dPhi -= TMath::TwoPi();
 457   if(dPhi<(-1.*TMath::Pi()))
 458     dPhi += TMath::TwoPi();
 459   Double_t dEta = jet1->Eta() - jet2->Eta();
 460   Double_t dR = TMath::Sqrt(dPhi*dPhi+dEta*dEta);
 461   return dR;
 462 }
 463
 464 //________________________________________________________________________
 465 Double_t AliAnalysisTaskEmcalJetTagger::GetDeltaPhi(const AliEmcalJet* jet1, const AliEmcalJet* jet2) {
 466   //
 467   // Calculate azimuthal angle between the axises of the jets
 468   //
 469
 470   return GetDeltaPhi(jet1->Phi(),jet2->Phi());
 471
 472 }
 473
 474 //________________________________________________________________________
 475 Double_t AliAnalysisTaskEmcalJetTagger::GetDeltaPhi(Double_t phi1,Double_t phi2) {
 476   //
 477   // Calculate azimuthal angle between the axises of the jets
 478   //
 479
 480   Double_t dPhi = phi1-phi2;
 481
 482   if(dPhi <-0.5*TMath::Pi())  dPhi += TMath::TwoPi();
 483   if(dPhi > 1.5*TMath::Pi())  dPhi -= TMath::TwoPi();
 484
 485   return dPhi;
 486 }
 487
 488 //________________________________________________________________________
 489 Double_t AliAnalysisTaskEmcalJetTagger::GetFractionSharedPt(const AliEmcalJet *jet1, const AliEmcalJet *jet2) const
 490 {
 491   //
 492   // Get fraction of shared pT between matched full and charged jet
 493   // Uses charged jet pT as baseline: fraction = \Sum_{const,full jet} pT,const,i / pT,jet,ch
 494   //
 495
 496   Double_t fraction = 0.;
 497   Double_t jetPt2 = jet2->Pt();
 498
 499   if(jetPt2>0) {
 500
 501     AliDebug(2,Form("%s: nConstituents_jet1: %d, nConstituents_jet2: %d, tracks_jet1: %d  tracks_jet2: %d clusters_jet1: %d clusters_jet2: %d",GetName(),jet1->GetNumberOfConstituents(),jet2->GetNumberOfConstituents(),jet1->GetNumberOfTracks(),jet2->GetNumberOfTracks(),jet1->GetNumberOfClusters(),jet2->GetNumberOfClusters()));
 502
 503     Double_t sumPt = 0.;
 504     AliVParticle *vpf = 0x0;
 505     Int_t iFound = 0;
 506     for(Int_t icc=0; icc<jet2->GetNumberOfTracks(); icc++) {
 507       Int_t idx = (Int_t)jet2->TrackAt(icc);
 508       iFound = 0;
 509       for(Int_t icf=0; icf<jet1->GetNumberOfTracks(); icf++) {
 510         if(idx == jet1->TrackAt(icf) && iFound==0 ) {
 511           iFound=1;
 512           vpf = static_cast<AliVParticle*>(jet1->TrackAt(icf, fTracks));
 513           if(vpf) sumPt += vpf->Pt();
 514           continue;
 515         }
 516       }
 517     }
 518     fraction = sumPt/jetPt2;
 519   } else fraction = -1.;
 520   return fraction;
 521 }
 522
 523 //________________________________________________________________________
 524 Bool_t AliAnalysisTaskEmcalJetTagger::RetrieveEventObjects() {
 525   //
 526   // retrieve event objects
 527   //
 528
 529   if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
 530     return kFALSE;
 531
 532   return kTRUE;
 533
 534 }
 535
 536 //_______________________________________________________________________
 537 void AliAnalysisTaskEmcalJetTagger::Terminate(Option_t *)
 538 {
 539   // Called once at the end of the analysis.
 540 }
 541