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 "AliAnalysisUtils.h"
  25 #include "AliEmcalParticle.h"
  26 #include "AliMCEvent.h"
  27 #include "AliGenPythiaEventHeader.h"
  28 #include "AliAODMCHeader.h"
  29 #include "AliMCEvent.h"
  30 #include "AliAnalysisManager.h"
  31 #include "AliJetContainer.h"
  32
  33 #include "AliAnalysisTaskEmcalDiJetBase.h"
  34
  35 ClassImp(AliAnalysisTaskEmcalDiJetBase)
  36
  37 //________________________________________________________________________
  38 AliAnalysisTaskEmcalDiJetBase::AliAnalysisTaskEmcalDiJetBase() :
  39   AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalDiJetBase", kTRUE),
  40   fDebug(kFALSE),
  41   fJetCorrelationType(kCorrelateAll),
  42   fJetFullChargedMatchingType(kFraction),
  43   fTriggerClass(""),
  44   fContainerCharged(1),
  45   fContainerFull(0),
  46   fContainerChargedMC(3),
  47   fContainerFullMC(2),
  48   fRhoType(0),
  49   fRhoChVal(0),
  50   fRhoFullVal(0),
  51   fDoChargedCharged(kTRUE),
  52   fDoFullCharged(kTRUE),
  53   fDoFullFull(kFALSE),
  54   fUseAnaUtils(kTRUE),
  55   fAnalysisUtils(0),
  56   fPtMinTriggerJet(10.),
  57   fMinFractionShared(0.5),
  58   fMatchingDone(kFALSE),
  59   faFullFracIndex(0),
  60   faFullFracIndexMC(0),
  61   fIsPythiaPtHard(0),
  62   fPythiaHeader(0),
  63   fPtHard(0),
  64   fPtHardBin(0),
  65   fNTrials(0),
  66   fhNEvents(0),
  67   fHistTrials(0),
  68   fHistTrialsSelEvents(0),
  69   fHistXsection(0),
  70   fHistEvents(0)
  71 {
  72   // Default constructor.
  73
  74
  75   SetMakeGeneralHistograms(kTRUE);
  76
  77 }
  78
  79 //________________________________________________________________________
  80 AliAnalysisTaskEmcalDiJetBase::AliAnalysisTaskEmcalDiJetBase(const char *name) :
  81   AliAnalysisTaskEmcalJet(name, kTRUE),
  82   fDebug(kFALSE),
  83   fJetCorrelationType(kCorrelateAll),
  84   fJetFullChargedMatchingType(kFraction),
  85   fTriggerClass(""),
  86   fContainerCharged(1),
  87   fContainerFull(0),
  88   fContainerChargedMC(3),
  89   fContainerFullMC(2),
  90   fRhoType(0),
  91   fRhoChVal(0),
  92   fRhoFullVal(0),
  93   fDoChargedCharged(kTRUE),
  94   fDoFullCharged(kTRUE),
  95   fDoFullFull(kFALSE),
  96   fUseAnaUtils(kTRUE),
  97   fAnalysisUtils(0),
  98   fPtMinTriggerJet(10.),
  99   fMinFractionShared(0.5),
 100   fMatchingDone(kFALSE),
 101   faFullFracIndex(0),
 102   faFullFracIndexMC(0),
 103   fIsPythiaPtHard(0),
 104   fPythiaHeader(0),
 105   fPtHard(0),
 106   fPtHardBin(0),
 107   fNTrials(0),
 108   fhNEvents(0),
 109   fHistTrials(0),
 110   fHistTrialsSelEvents(0),
 111   fHistXsection(0),
 112   fHistEvents(0)
 113 {
 114   // Standard constructor.
 115
 116   SetMakeGeneralHistograms(kTRUE);
 117
 118 }
 119
 120 //________________________________________________________________________
 121 AliAnalysisTaskEmcalDiJetBase::~AliAnalysisTaskEmcalDiJetBase()
 122 {
 123   // Destructor.
 124 }
 125
 126 //----------------------------------------------------------------------------------------------
 127 void AliAnalysisTaskEmcalDiJetBase::InitOnce() {
 128   //
 129   // only initialize once
 130   //
 131
 132   // Initialize analysis util class for vertex selection
 133   if(fUseAnaUtils) {
 134     fAnalysisUtils = new AliAnalysisUtils();
 135     fAnalysisUtils->SetMinVtxContr(2);
 136     fAnalysisUtils->SetMaxVtxZ(10.);
 137   }
 138 }
 139
 140 //________________________________________________________________________
 141 Bool_t AliAnalysisTaskEmcalDiJetBase::SelectEvent() {
 142   //
 143   // Decide if event should be selected for analysis
 144   //
 145
 146   fhNEvents->Fill(0.5);
 147
 148   if(!fAnalysisUtils && fUseAnaUtils)
 149     InitOnce();
 150
 151
 152  if(fAnalysisUtils) {
 153     if(!fAnalysisUtils->IsVertexSelected2013pA(InputEvent()))
 154       return kFALSE;
 155
 156     fhNEvents->Fill(2.5);
 157
 158     if(fAnalysisUtils->IsPileUpEvent(InputEvent()))
 159       return kFALSE;
 160   }
 161   else{
 162     if(fUseAnaUtils)
 163       AliError("fAnalysisUtils not initialized. Call AliAnalysisTaskEmcalJetTriggerQA::InitOnce()");
 164   }
 165
 166   fhNEvents->Fill(3.5);
 167
 168   if(!fTriggerClass.IsNull()) {
 169     //Check if requested trigger was fired
 170     TString firedTrigClass = InputEvent()->GetFiredTriggerClasses();
 171     if(!firedTrigClass.Contains(fTriggerClass))
 172       return kFALSE;
 173   }
 174
 175   fhNEvents->Fill(1.5);
 176
 177   return kTRUE;
 178
 179 }
 180
 181 //________________________________________________________________________
 182 void AliAnalysisTaskEmcalDiJetBase::UserCreateOutputObjects()
 183 {
 184   // Create user output.
 185
 186   InitOnce();
 187
 188   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 189
 190   Bool_t oldStatus = TH1::AddDirectoryStatus();
 191   TH1::AddDirectory(kFALSE);
 192
 193   fhNEvents = new TH1F("fhNEvents","fhNEvents;selection;N_{evt}",5,0,5);
 194   fOutput->Add(fhNEvents);
 195
 196
 197   //Pythia info
 198
 199   fHistTrials = new TH1F("fHistTrials", "fHistTrials", 12, 0, 12);
 200   fHistTrials->GetXaxis()->SetTitle("p_{T} hard bin");
 201   fHistTrials->GetYaxis()->SetTitle("trials");
 202   fOutput->Add(fHistTrials);
 203
 204   fHistTrialsSelEvents = new TH1F("fHistTrialsSelEvents", "fHistTrialsSelEvents", 12, 0, 12);
 205   fHistTrialsSelEvents->GetXaxis()->SetTitle("p_{T} hard bin");
 206   fHistTrialsSelEvents->GetYaxis()->SetTitle("trials");
 207   fOutput->Add(fHistTrialsSelEvents);
 208
 209   fHistXsection = new TProfile("fHistXsection", "fHistXsection", 12, 0, 12);
 210   fHistXsection->GetXaxis()->SetTitle("p_{T} hard bin");
 211   fHistXsection->GetYaxis()->SetTitle("xsection");
 212   fOutput->Add(fHistXsection);
 213
 214   fHistEvents = new TH1F("fHistEvents", "fHistEvents", 12, 0, 12);
 215   fHistEvents->GetXaxis()->SetTitle("p_{T} hard bin");
 216   fHistEvents->GetYaxis()->SetTitle("total events");
 217   fOutput->Add(fHistEvents);
 218
 219
 220   TH1::AddDirectory(oldStatus);
 221
 222   PostData(1, fOutput); // Post data for ALL output slots > 0 here.
 223 }
 224
 225 //________________________________________________________________________
 226 Bool_t AliAnalysisTaskEmcalDiJetBase::IsSameJet(const Int_t ijt, const Int_t ija, const Int_t type, const Bool_t isMC) {
 227    //check if two jets are the same one
 228
 229    Bool_t bSame = kFALSE;
 230
 231    if(type<2 && ijt==ija)
 232      bSame = kTRUE;
 233    if(type==2) {
 234
 235      if(fJetFullChargedMatchingType==kFraction) {
 236        if(isMC && faFullFracIndexMC.At(ijt)==ija)
 237          bSame = kTRUE;
 238        else if(!isMC && faFullFracIndex.At(ijt)==ija)
 239          bSame = kTRUE;
 240      }
 241      else if(fJetFullChargedMatchingType==kGeo) {
 242        Int_t contFull = fContainerFull;
 243        Int_t contChar = fContainerCharged;
 244        if(isMC) {
 245          contFull = fContainerFullMC;
 246          contChar = fContainerChargedMC;
 247        }
 248        AliEmcalJet *fullJet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ijt, contFull));
 249        AliEmcalJet *chJet   = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ija, contChar));
 250        AliEmcalJet *matchJet = fullJet->ClosestJet();
 251        if(chJet==matchJet)
 252          bSame = kTRUE;
 253      }
 254      else if(fJetFullChargedMatchingType==kNoMatching) {
 255        return kFALSE;
 256      }
 257      else{
 258        AliWarning(Form("%s: matching type unknown", GetName()));
 259      }
 260
 261
 262    }
 263
 264    return bSame;
 265  }
 266
 267
 268 //________________________________________________________________________
 269 Double_t AliAnalysisTaskEmcalDiJetBase::GetJetPt(const AliEmcalJet *jet, const Int_t type) {
 270
 271   if(!jet) return -99;
 272
 273   if(type==0)
 274     return jet->Pt() - fRhoFullVal * jet->Area();
 275   else if(type==1)
 276     return jet->Pt() - fRhoChVal * jet->Area();
 277   else
 278     return jet->Pt();
 279
 280 }
 281
 282 //________________________________________________________________________
 283 Double_t AliAnalysisTaskEmcalDiJetBase::GetZ(const AliVParticle *trk, const AliEmcalJet *jet)          const
 284 {
 285   // Get Z of constituent trk
 286
 287   return GetZ(trk->Px(),trk->Py(),trk->Pz(),jet->Px(),jet->Py(),jet->Pz());
 288 }
 289
 290 //________________________________________________________________________
 291 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
 292 {
 293   //
 294   // Get the z of a constituent inside of a jet
 295   //
 296
 297   Double_t pJetSq = jetPx*jetPx+jetPy*jetPy+jetPz*jetPz;
 298
 299   if(pJetSq>0.)
 300     return (trkPx*jetPx+trkPy*jetPy+trkPz*jetPz)/pJetSq;
 301   else {
 302     AliWarning(Form("%s: strange, pjet*pjet seems to be zero pJetSq: %f",GetName(), pJetSq));
 303     return 0;
 304   }
 305 }
 306
 307 //________________________________________________________________________
 308 Double_t AliAnalysisTaskEmcalDiJetBase::GetDeltaR(const AliEmcalJet* jet1, const AliEmcalJet* jet2) const {
 309   //
 310   // Helper function to calculate the distance between two jets
 311   //
 312
 313   Double_t dPhi = jet1->Phi() - jet2->Phi();
 314   if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 315   if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 316   Double_t dEta = jet1->Eta() - jet2->Eta();
 317   Double_t dR = TMath::Sqrt(dPhi*dPhi+dEta*dEta);
 318   return dR;
 319 }
 320
 321 //________________________________________________________________________
 322 Double_t AliAnalysisTaskEmcalDiJetBase::GetDeltaPhi(const AliEmcalJet* jet1, const AliEmcalJet* jet2) {
 323   //
 324   // Calculate azimuthal angle between the axises of the jets
 325   //
 326
 327   return GetDeltaPhi(jet1->Phi(),jet2->Phi());
 328
 329 }
 330
 331 //________________________________________________________________________
 332 Double_t AliAnalysisTaskEmcalDiJetBase::GetDeltaPhi(Double_t phi1,Double_t phi2) {
 333   //
 334   // Calculate azimuthal angle between the axises of the jets
 335   //
 336
 337   Double_t dPhi = phi1-phi2;
 338
 339   if(dPhi > 2*TMath::Pi())    dPhi -= TMath::TwoPi();
 340   if(dPhi <-2*TMath::Pi())    dPhi += TMath::TwoPi();
 341   if(dPhi <-0.5*TMath::Pi())  dPhi += TMath::TwoPi();
 342   if(dPhi > 1.5*TMath::Pi())  dPhi -= TMath::TwoPi();
 343
 344   return dPhi;
 345 }
 346
 347 //________________________________________________________________________
 348 Double_t AliAnalysisTaskEmcalDiJetBase::GetFractionSharedPt(const AliEmcalJet *jetFull, const AliEmcalJet *jetCharged) const
 349 {
 350   //
 351   // Get fraction of shared pT between matched full and charged jet
 352   // Uses charged jet pT as baseline: fraction = \Sum_{const,full jet} pT,const,i / pT,jet,ch
 353   //
 354
 355   Double_t fraction = 0.;
 356   Double_t jetPtCh = jetCharged->Pt();
 357
 358   if(jetPtCh>0) {
 359
 360     if(fDebug>10)  AliInfo(Form("%s: nConstituents: %d, ch: %d  chne: %d ne: %d",GetName(),jetFull->GetNumberOfConstituents(),jetCharged->GetNumberOfTracks(),jetFull->GetNumberOfTracks(),jetFull->GetNumberOfClusters()));
 361
 362     Double_t sumPt = 0.;
 363     AliVParticle *vpf = 0x0;
 364     Int_t iFound = 0;
 365     for(Int_t icc=0; icc<jetCharged->GetNumberOfTracks(); icc++) {
 366       Int_t idx = (Int_t)jetCharged->TrackAt(icc);
 367       iFound = 0;
 368       for(Int_t icf=0; icf<jetFull->GetNumberOfTracks(); icf++) {
 369         if(idx == jetFull->TrackAt(icf) && iFound==0 ) {
 370           iFound=1;
 371           vpf = static_cast<AliVParticle*>(jetFull->TrackAt(icf, fTracks));
 372           if(!vpf) continue;
 373           if(vpf->Charge()!=0)
 374             sumPt += vpf->Pt();
 375           continue;
 376         }
 377       }
 378     }
 379
 380     fraction = sumPt/jetPtCh;
 381   }
 382
 383   if(fDebug>10) AliInfo(Form("%s: charged shared fraction: %.2f",GetName(),fraction));
 384
 385   return fraction;
 386
 387 }
 388
 389 //________________________________________________________________________
 390 void AliAnalysisTaskEmcalDiJetBase::MatchJetsGeo(Int_t cFull, Int_t cCharged,
 391                                                  Int_t iDebug, Float_t maxDist, Int_t type) {
 392
 393   //
 394   // Match the full jets to the corresponding charged jets
 395   // Translation of AliAnalysisHelperJetTasks::GetClosestJets to AliEmcalJet objects
 396   // 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)
 397
 398   const int kMode = 3;
 399
 400   const Int_t nChJets   = GetNJets(cCharged);
 401   const Int_t nFullJets = GetNJets(cFull);
 402
 403   if(nFullJets==0 || nChJets==0) return;
 404
 405   TArrayI faChNeMatchIndex;
 406   faChNeMatchIndex.Set(nChJets+1);
 407
 408   TArrayI faChMatchIndex;
 409   faChMatchIndex.Set(nFullJets+1);
 410
 411   static TArrayS iFlag(nChJets*nFullJets);
 412   if(iFlag.GetSize()<(nChJets*nFullJets)){
 413     iFlag.Set(nChJets*nFullJets+1);
 414   }
 415   iFlag.Reset(0);
 416
 417   AliJetContainer *contCh = GetJetContainer(cCharged);
 418
 419   // find the closest distance to the full jet
 420   for(int ifu = 0;ifu<nFullJets;ifu++){
 421
 422     AliEmcalJet *fullJet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ifu, cFull));
 423     if(!fullJet) continue;
 424
 425     Float_t dist = maxDist;
 426
 427     for(int ich = 0;ich<nChJets;ich++){
 428       AliEmcalJet *chJet = 0x0;
 429       if(type==0)
 430         chJet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ich, cCharged));
 431       else {
 432         chJet = static_cast<AliEmcalJet*>(GetJetFromArray(ich, cCharged));
 433         if(!chJet) continue;;
 434         if(type>0) {
 435           if(chJet->Eta()<(contCh->GetJetEtaMin()-0.1) || chJet->Eta()>(contCh->GetJetEtaMax()+0.1))
 436             continue;
 437           if(type==2) {
 438             if(chJet->Phi()<(contCh->GetJetPhiMin()-0.1) || chJet->Phi()>(contCh->GetJetPhiMax()+0.1))
 439               continue;
 440           }
 441         }
 442       }
 443       if(!chJet)
 444         continue;
 445
 446       Double_t frac = GetFractionSharedPt(fullJet,chJet);
 447       Double_t dR = GetDeltaR(fullJet,chJet);
 448       if(dR<dist){
 449         faChMatchIndex[ifu]=ich;
 450         dist = dR;
 451       }
 452       if(iDebug>10) Printf("Distance (%d)--(%d) %3.3f frac:%.2f",ifu,ich,dR,frac);
 453     }
 454     if(faChMatchIndex[ifu]>=0) iFlag[ifu*nChJets+faChMatchIndex[ifu]]+=1;//ich closest to ifu
 455     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]]);
 456
 457     // reset...
 458     faChMatchIndex[ifu]=-1;
 459
 460
 461   }//full jet loop
 462
 463
 464   // other way around
 465   for(int ich = 0;ich<nChJets;ich++){
 466     AliEmcalJet *chJet = 0x0;
 467     if(type==0)
 468       chJet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ich, cCharged));
 469     else {
 470       chJet = static_cast<AliEmcalJet*>(GetJetFromArray(ich, cCharged));
 471       if(!chJet) continue;;
 472       if(type>0) {
 473         if(chJet->Eta()<(contCh->GetJetEtaMin()-0.1) || chJet->Eta()>(contCh->GetJetEtaMax()+0.1))
 474           continue;
 475         if(type==2) {
 476           if(chJet->Phi()<(contCh->GetJetPhiMin()-0.1) || chJet->Phi()>(contCh->GetJetPhiMax()+0.1))
 477             continue;
 478         }
 479       }
 480     }
 481     if(!chJet)
 482       continue;
 483
 484     Float_t dist = maxDist;
 485     for(int ifu = 0;ifu<nFullJets;ifu++){
 486       AliEmcalJet *fullJet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ifu, cFull));
 487       if(!fullJet)
 488         continue;
 489
 490       Double_t dR = GetDeltaR(fullJet,chJet);
 491       if(dR<dist){
 492         faChNeMatchIndex[ich]=ifu;
 493         dist = dR;
 494       }
 495     }
 496     if(faChNeMatchIndex[ich]>=0) iFlag[faChNeMatchIndex[ich]*nChJets+ich]+=2;//ifu closest to ich
 497     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]);
 498
 499     // reset...
 500     faChNeMatchIndex[ich]=-1;
 501
 502   }
 503
 504
 505   // check for "true" correlations
 506   for(int ifu = 0;ifu<nFullJets;ifu++){
 507     for(int ich = 0;ich<nChJets;ich++){
 508       if(iDebug>10) AliInfo(Form("%s: Flag[%d][%d] %d ",GetName(),ifu,ich,iFlag[ifu*nChJets+ich]));
 509
 510       if(kMode==3){
 511         // we have a uniqe correlation
 512         if(iFlag[ifu*nChJets+ich]==3){
 513
 514           AliEmcalJet *chJet = static_cast<AliEmcalJet*>(GetJetFromArray(ich, cCharged));
 515           AliEmcalJet *fullJet = static_cast<AliEmcalJet*>(GetJetFromArray(ifu, cFull));
 516           Double_t dR = GetDeltaR(fullJet,chJet);
 517
 518           if(iDebug>10) Printf("closest jets %d  %d  dR =  %f",ich,ifu,dR);
 519
 520           chJet->SetClosestJet(fullJet,dR);
 521           fullJet->SetClosestJet(chJet,dR);
 522
 523         }
 524       }
 525     }
 526   }
 527
 528   fMatchingDone = kTRUE;
 529
 530 }
 531
 532 //________________________________________________________________________
 533 void AliAnalysisTaskEmcalDiJetBase::SetChargedFractionIndex() {
 534
 535   // take each full jet and set the index of the charged jet with the largest shared charged fraction
 536
 537   const Int_t nJetsCh   = GetNJets(fContainerCharged);
 538   const Int_t nJetsFull = GetNJets(fContainerFull);
 539   faFullFracIndex.Set(nJetsFull+1);
 540   faFullFracIndex.Reset(-1);
 541
 542   AliJetContainer  *cont = GetJetContainer(fContainerFull);
 543   Float_t radius =  cont->GetJetRadius();
 544
 545   for(Int_t ifu = 0; ifu<nJetsFull; ifu++) {
 546     Double_t frac = 0.;
 547     Double_t dist = 10.;
 548     AliEmcalJet *jetFull = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ifu, fContainerFull));
 549     if(!jetFull) {
 550       faFullFracIndex.SetAt(-1,ifu);
 551       continue;
 552     }
 553
 554     for(Int_t ich = 0; ich<nJetsCh; ich++) {
 555       AliEmcalJet *jetCh = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ich, fContainerCharged));
 556       if(!jetCh)
 557         continue;
 558       Double_t tmpFrac = GetFractionSharedPt(jetFull,jetCh);
 559       dist = GetDeltaR(jetFull,jetCh);
 560       if(tmpFrac>frac && dist<radius) {
 561         frac = tmpFrac;
 562         faFullFracIndex.SetAt(ich,ifu);
 563       }
 564     }
 565
 566   }
 567
 568 }
 569
 570 //________________________________________________________________________
 571 void AliAnalysisTaskEmcalDiJetBase::SetChargedFractionIndexMC() {
 572
 573   // take each full jet and set the index of the charged jet with the largest shared charged fraction
 574
 575   const Int_t nJetsCh   = GetNJets(fContainerChargedMC);
 576   const Int_t nJetsFull = GetNJets(fContainerFullMC);
 577   faFullFracIndexMC.Set(nJetsFull);
 578   faFullFracIndexMC.Reset(-1);
 579
 580   AliJetContainer  *cont = GetJetContainer(fContainerFullMC);
 581   Float_t radius =  cont->GetJetRadius();
 582
 583   for(Int_t ifu = 0; ifu<nJetsFull; ifu++) {
 584     Double_t frac = 0.;
 585     Double_t dist = 10.;
 586     AliEmcalJet *jetFull = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ifu, fContainerFullMC));
 587     if(!jetFull) {
 588       faFullFracIndexMC.SetAt(-1,ifu);
 589       continue;
 590     }
 591
 592     for(Int_t ich = 0; ich<nJetsCh; ich++) {
 593       AliEmcalJet *jetCh = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ich, fContainerChargedMC));
 594       if(!jetCh)
 595         continue;
 596       Double_t tmpFrac = GetFractionSharedPt(jetFull,jetCh);
 597       dist = GetDeltaR(jetFull,jetCh);
 598       if(tmpFrac>frac && dist<radius) {
 599         frac = tmpFrac;
 600         faFullFracIndexMC.SetAt(ich,ifu);
 601       }
 602     }
 603
 604   }
 605
 606 }
 607
 608 //_______________________________________________________________________
 609 AliEmcalJet* AliAnalysisTaskEmcalDiJetBase::GetLeadingJetOppositeHemisphere(const Int_t type, const Int_t typea, const AliEmcalJet *jetTrig) {
 610
 611   // Get leading jet in opposite hemisphere from trigger jet
 612   // type = correlation type
 613   // typea = container of associated jets
 614
 615   Int_t nJetsAssoc = GetNJets(typea);
 616   Double_t ptLead = -999;
 617   Int_t    iJetLead = -1;
 618   for(Int_t ija=0; ija<nJetsAssoc; ija++) {
 619
 620     AliEmcalJet *jetAssoc = NULL;
 621     if(type==0) {
 622       jetAssoc = static_cast<AliEmcalJet*>(GetJetFromArray(ija, typea));
 623       if(TMath::Abs(jetAssoc->Eta())>0.5)
 624         jetAssoc = NULL;
 625     }
 626     else
 627       jetAssoc = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ija, typea));
 628
 629     if(!jetAssoc)
 630       continue;
 631
 632     Double_t dPhi = GetDeltaPhi(jetTrig,jetAssoc);
 633     Double_t phiMin = 0.5*TMath::Pi();
 634     Double_t phiMax = 1.5*TMath::Pi();
 635     if(dPhi<phiMin || dPhi>phiMax)
 636       continue;
 637
 638     Double_t jetAssocPt = GetJetPt(jetAssoc,typea);
 639
 640     if(jetAssocPt>ptLead) {
 641       ptLead = jetAssocPt;
 642       iJetLead = ija;
 643     }
 644
 645   }
 646
 647   AliEmcalJet *jetAssocLead = static_cast<AliEmcalJet*>(GetJetFromArray(iJetLead, typea));
 648
 649   return jetAssocLead;
 650
 651 }
 652
 653 //________________________________________________________________________
 654 Bool_t AliAnalysisTaskEmcalDiJetBase::PythiaInfoFromFile(const char* currFile, Float_t &xsec, Float_t &trials, Int_t &pthard)
 655 {
 656   //
 657   // Get the cross section and the trails either from pyxsec.root or from pysec_hists.root
 658   // Get the pt hard bin from the file path
 659   // This is to called in Notify and should provide the path to the AOD/ESD file
 660   // (Partially copied from AliAnalysisHelperJetTasks)
 661
 662   TString file(currFile);
 663   xsec = 0;
 664   trials = 1;
 665
 666   if(file.Contains(".zip#")){
 667     Ssiz_t pos1 = file.Index("root_archive",12,0,TString::kExact);
 668     Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
 669     Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
 670     file.Replace(pos+1,pos2-pos1,"");
 671   }
 672   else {
 673     // not an archive take the basename....
 674     file.ReplaceAll(gSystem->BaseName(file.Data()),"");
 675   }
 676   //Printf("%s",file.Data());
 677
 678   // Get the pt hard bin
 679   TString strPthard(file);
 680
 681   strPthard.Remove(strPthard.Last('/'));
 682   strPthard.Remove(strPthard.Last('/'));
 683   if (strPthard.Contains("AOD")) strPthard.Remove(strPthard.Last('/'));
 684   strPthard.Remove(0,strPthard.Last('/')+1);
 685   if (strPthard.IsDec())
 686     pthard = strPthard.Atoi();
 687   else
 688     AliWarning(Form("Could not extract file number from path %s", strPthard.Data()));
 689
 690   TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")); // problem that we cannot really test the existance of a file in a archive so we have to lvie with open error message from root
 691   if(!fxsec){
 692     // next trial fetch the histgram file
 693     fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
 694     if(!fxsec){
 695         // not a severe condition but inciate that we have no information
 696       return kFALSE;
 697     }
 698     else{
 699       // find the tlist we want to be independtent of the name so use the Tkey
 700       TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
 701       if(!key){
 702         fxsec->Close();
 703         return kFALSE;
 704       }
 705       TList *list = dynamic_cast<TList*>(key->ReadObj());
 706       if(!list){
 707         fxsec->Close();
 708         return kFALSE;
 709       }
 710       xsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
 711       trials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
 712       fxsec->Close();
 713     }
 714   } // no tree pyxsec.root
 715   else {
 716     TTree *xtree = (TTree*)fxsec->Get("Xsection");
 717     if(!xtree){
 718       fxsec->Close();
 719       return kFALSE;
 720     }
 721     UInt_t   ntrials  = 0;
 722     Double_t  xsection  = 0;
 723     xtree->SetBranchAddress("xsection",&xsection);
 724     xtree->SetBranchAddress("ntrials",&ntrials);
 725     xtree->GetEntry(0);
 726     trials = ntrials;
 727     xsec = xsection;
 728     fxsec->Close();
 729   }
 730   return kTRUE;
 731 }
 732
 733 //________________________________________________________________________
 734 Bool_t AliAnalysisTaskEmcalDiJetBase::UserNotify()
 735 {
 736   // Get pythia info
 737
 738   if (!fIsPythiaPtHard) {
 739
 740     return kTRUE;
 741   }
 742
 743   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 744   if (!tree) {
 745     AliError(Form("%s - UserNotify: No current tree!",GetName()));
 746     return kFALSE;
 747   }
 748
 749   Float_t xsection = 0;
 750   Float_t trials   = 0;
 751   Int_t   pthard   = 0;
 752
 753   TFile *curfile = tree->GetCurrentFile();
 754   if (!curfile) {
 755     AliError(Form("%s - UserNotify: No current file!",GetName()));
 756     return kFALSE;
 757   }
 758
 759   TChain *chain = dynamic_cast<TChain*>(tree);
 760   if (chain)
 761     tree = chain->GetTree();
 762
 763   Int_t nevents = tree->GetEntriesFast();
 764
 765   PythiaInfoFromFile(curfile->GetName(), xsection, trials, pthard);
 766
 767   fPtHardBin = pthard;
 768   // fXsec = xsection;
 769   // fTrials = trials;
 770
 771   fHistTrials->Fill(pthard, trials);
 772   fHistXsection->Fill(pthard, xsection);
 773   fHistEvents->Fill(pthard, nevents);
 774
 775   return kTRUE;
 776 }
 777
 778 //________________________________________________________________________
 779 Bool_t AliAnalysisTaskEmcalDiJetBase::RetrieveEventObjects() {
 780   //
 781   // get charged jets
 782   //
 783
 784   if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
 785     return kFALSE;
 786
 787   if(fRhoType==0) {
 788     fRhoFullVal = 0.;
 789     fRhoChVal = 0.;
 790   }
 791   if(fRhoType==1) {
 792     fRhoFullVal = GetRhoVal(fContainerFull);
 793     fRhoChVal = GetRhoVal(fContainerCharged);
 794   }
 795
 796   if (MCEvent()) {
 797     fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
 798     if (!fPythiaHeader) {
 799       // Check if AOD
 800       AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
 801
 802       if (aodMCH) {
 803         for(UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++) {
 804           fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
 805           if (fPythiaHeader) break;
 806         }
 807       }
 808     }
 809   }
 810
 811   if (fPythiaHeader) {
 812     fPtHard = fPythiaHeader->GetPtHard();
 813
 814     fNTrials = fPythiaHeader->Trials();
 815
 816   }
 817
 818   return kTRUE;
 819
 820 }
 821
 822 //_______________________________________________________________________
 823 void AliAnalysisTaskEmcalDiJetBase::Terminate(Option_t *)
 824 {
 825   // Called once at the end of the analysis.
 826 }