PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalDiJetAna.cxx

   1 //
   2 // Dijet 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 "AliCentrality.h"
  32
  33 #include "AliAnalysisTaskEmcalDiJetAna.h"
  34
  35 ClassImp(AliAnalysisTaskEmcalDiJetAna)
  36
  37 //________________________________________________________________________
  38 AliAnalysisTaskEmcalDiJetAna::AliAnalysisTaskEmcalDiJetAna() :
  39   AliAnalysisTaskEmcalDiJetBase("AliAnalysisTaskEmcalDiJetAna"),
  40   fDoMatchFullCharged(kTRUE),
  41   fNKtBins(30),
  42   fNDiJetEtaBins(1),
  43   fNAjBins(1),
  44   fh2CentRhoCh(0),
  45   fh2CentRhoScaled(0),
  46   fh3PtEtaPhiJetFull(0),
  47   fh3PtEtaPhiJetCharged(0),
  48   fhnDiJetVarsFull(0),
  49   fhnDiJetVarsCh(0),
  50   fhnDiJetVarsFullCharged(0),
  51   fhnMatchingFullCharged(0),
  52   fh3PtTrigKt1Kt2Ch(0),
  53   fh3PtTrigKt1Kt2FuCh(0),
  54   fh3PtTrigDPhi1DPhi2Ch(0),
  55   fh3PtTrigDPhi1DPhi2FuCh(0)
  56 {
  57   // Default constructor.
  58
  59   for(Int_t i=0; i<4; i++) {
  60     fh3DiJetKtNEFPtAssoc[i]          = 0;
  61     fCentCorrPtAssocCh[i]            = 0;
  62     fCentCorrPtAssocFuCh[i]          = 0;
  63     fAjPtAssocCentCh[i]              = 0;
  64     fAjPtAssocCentFuCh[i]            = 0;
  65     fh3PtAssoc1PtAssoc2DPhi23Ch[i]   = 0;
  66     fh3PtAssoc1PtAssoc2DPhi23FuCh[i] = 0;
  67   }
  68
  69   SetMakeGeneralHistograms(kTRUE);
  70 }
  71
  72 //________________________________________________________________________
  73 AliAnalysisTaskEmcalDiJetAna::AliAnalysisTaskEmcalDiJetAna(const char *name) :
  74   AliAnalysisTaskEmcalDiJetBase(name),
  75   fDoMatchFullCharged(kTRUE),
  76   fNKtBins(30),
  77   fNDiJetEtaBins(1),
  78   fNAjBins(1),
  79   fh2CentRhoCh(0),
  80   fh2CentRhoScaled(0),
  81   fh3PtEtaPhiJetFull(0),
  82   fh3PtEtaPhiJetCharged(0),
  83   fhnDiJetVarsFull(0),
  84   fhnDiJetVarsCh(0),
  85   fhnDiJetVarsFullCharged(0),
  86   fhnMatchingFullCharged(0),
  87   fh3PtTrigKt1Kt2Ch(0),
  88   fh3PtTrigKt1Kt2FuCh(0),
  89   fh3PtTrigDPhi1DPhi2Ch(0),
  90   fh3PtTrigDPhi1DPhi2FuCh(0)
  91 {
  92   // Standard constructor.
  93
  94   for(Int_t i=0; i<4; i++) {
  95     fh3DiJetKtNEFPtAssoc[i]      = 0;
  96     fCentCorrPtAssocCh[i]        = 0;
  97     fCentCorrPtAssocFuCh[i]      = 0;
  98     fAjPtAssocCentCh[i]          = 0;
  99     fAjPtAssocCentFuCh[i]        = 0;
 100     fh3PtAssoc1PtAssoc2DPhi23Ch[i]   = 0;
 101     fh3PtAssoc1PtAssoc2DPhi23FuCh[i] = 0;
 102   }
 103
 104   SetMakeGeneralHistograms(kTRUE);
 105 }
 106
 107 //________________________________________________________________________
 108 AliAnalysisTaskEmcalDiJetAna::~AliAnalysisTaskEmcalDiJetAna()
 109 {
 110   // Destructor.
 111 }
 112
 113 //________________________________________________________________________
 114 Bool_t AliAnalysisTaskEmcalDiJetAna::RetrieveEventObjects() {
 115   //
 116   // retrieve event objects
 117   //
 118
 119   if (!AliAnalysisTaskEmcalDiJetBase::RetrieveEventObjects())
 120     return kFALSE;
 121
 122   return kTRUE;
 123 }
 124
 125 //________________________________________________________________________
 126 void AliAnalysisTaskEmcalDiJetAna::UserCreateOutputObjects()
 127 {
 128   // Create user output.
 129
 130   AliAnalysisTaskEmcalDiJetBase::UserCreateOutputObjects();
 131
 132   Bool_t oldStatus = TH1::AddDirectoryStatus();
 133   TH1::AddDirectory(kFALSE);
 134
 135   const Int_t nBinsCent = 100;
 136   Double_t minCent = 0.;
 137   Double_t maxCent = 100.;
 138
 139   const Int_t nBinsRho = 200;
 140   Double_t minRho = 0.;
 141   Double_t maxRho = 20.;
 142   fh2CentRhoCh = new TH2F("fh2CentRhoCh","fh2CentRhoCh;centrality;#rho_{ch}",nBinsCent,minCent,maxCent,nBinsRho,minRho,maxRho);
 143   fOutput->Add(fh2CentRhoCh);
 144   fh2CentRhoScaled = new TH2F("fh2CentRhoScaled","fh2CentRhoScaled;centrality;s_{EMC}#rho_{ch}",nBinsCent,minCent,maxCent,nBinsRho,minRho,maxRho);
 145   fOutput->Add(fh2CentRhoScaled);
 146
 147   const Int_t nBinsPt = 150;
 148   Double_t minPt = -20.;
 149   Double_t maxPt = 130.;
 150   const Int_t nBinsEta = 40;
 151   Double_t minEta = -1.;
 152   Double_t maxEta = 1.;
 153   const Int_t nBinsPhi = 18*6;
 154   Double_t minPhi = 0.;
 155   Double_t maxPhi = TMath::TwoPi();
 156
 157   fh3PtEtaPhiJetFull = new TH3F("fh3PtEtaPhiJetFull","fh3PtEtaPhiJetFull;#it{p}_{T}^{jet};#eta;#varphi",nBinsPt,minPt,maxPt,nBinsEta,minEta,maxEta,nBinsPhi,minPhi,maxPhi);
 158   fOutput->Add(fh3PtEtaPhiJetFull);
 159
 160   fh3PtEtaPhiJetCharged = new TH3F("fh3PtEtaPhiJetCharged","fh3PtEtaPhiJetCharged;#it{p}_{T}^{jet};#eta;#varphi",nBinsPt,minPt,maxPt,nBinsEta,minEta,maxEta,nBinsPhi,minPhi,maxPhi);
 161   fOutput->Add(fh3PtEtaPhiJetCharged);
 162
 163   const Int_t nBinsSparse0 = 7;
 164   const Int_t nBinsPtW      = 30;
 165   const Int_t nBinsDPhi     = 72;
 166   const Int_t nBinsKt       = fNKtBins;
 167   const Int_t nBinsDiJetEta = fNDiJetEtaBins;//40;
 168   const Int_t nBinsCentr    = fNcentBins;
 169   const Int_t nBinsAj       = fNAjBins;//20
 170   const Int_t nBins0[nBinsSparse0] = {nBinsPtW,nBinsPtW,nBinsDPhi,nBinsKt,nBinsDiJetEta,nBinsCentr,nBinsAj};
 171   //pT1, pT2, deltaPhi, kT
 172   const Double_t xmin0[nBinsSparse0]  = {  minPt, minPt, -0.5*TMath::Pi(),   0.,-1.,0.  , 0.};
 173   const Double_t xmax0[nBinsSparse0]  = {  maxPt, maxPt,  1.5*TMath::Pi(), 120., 1.,100., 1.};
 174   const Double_t centArrayBins[8] = {0.,2.,5.,10.,20.,40.,60.,100.};
 175
 176   if(fDoChargedCharged) {
 177     fhnDiJetVarsCh = new THnSparseF("fhnDiJetVarsCh",
 178                                     "fhnDiJetVarsCh;#it{p}_{T,1} (GeV/#it{c});#it{p}_{T,2} (GeV/#it{c});#Delta#varphi;#it{k}_{T} = #it{p}_{T,1}sin(#Delta#varphi) (GeV/#it{c});(#eta_{1}+#eta_{2})/2);centrality;#it{A}_{j}",
 179                                     nBinsSparse0,nBins0,xmin0,xmax0);
 180     if(fNcentBins==7) fhnDiJetVarsCh->SetBinEdges(5,centArrayBins);
 181     fOutput->Add(fhnDiJetVarsCh);
 182   }
 183
 184   if(fDoFullCharged) {
 185     fhnDiJetVarsFullCharged = new THnSparseF("fhnDiJetVarsFullCharged",
 186                                 "fhnDiJetVarsFullCharged;#it{p}_{T,1} (GeV/#it{c});#it{p}_{T,2} (GeV/#it{c});#Delta#varphi;#it{k}_{T} = #it{p}_{T,1}sin(#Delta#varphi) (GeV/#it{c});(#eta_{1}+#eta_{2})/2);centrality;#it{A}_{j}",
 187                                 nBinsSparse0,nBins0,xmin0,xmax0);
 188     if(fNcentBins==7) fhnDiJetVarsFullCharged->SetBinEdges(5,centArrayBins);
 189     fOutput->Add(fhnDiJetVarsFullCharged);
 190   }
 191
 192   if(fDoFullFull) {
 193     fhnDiJetVarsFull = new THnSparseF("fhnDiJetVarsFull",
 194                                     "fhnDiJetVarsFull;#it{p}_{T,1} (GeV/#it{c});#it{p}_{T,2} (GeV/#it{c});#Delta#varphi;#it{k}_{T} = #it{p}_{T,1}sin(#Delta#varphi) (GeV/#it{c});(#eta_{1}+#eta_{2})/2);centrality;#it{A}_{j}",
 195                                     nBinsSparse0,nBins0,xmin0,xmax0);
 196     fOutput->Add(fhnDiJetVarsFull);
 197   }
 198
 199   fh3PtTrigKt1Kt2Ch = new TH3F("fh3PtTrigKt1Kt2Ch","fh3PtTrigKt1Kt2Ch;#it{p}_{T,1} (GeV/#it{c});#it{k}_{T,1};#it{k}_{T,2}",nBinsPt,minPt,maxPt,nBinsKt,0.,120.,nBinsKt,0.,120.);
 200   fOutput->Add(fh3PtTrigKt1Kt2Ch);
 201
 202   fh3PtTrigKt1Kt2FuCh = new TH3F("fh3PtTrigKt1Kt2FuCh","fh3PtTrigKt1Kt2FuCh;#it{p}_{T,1} (GeV/#it{c});#it{k}_{T,1};#it{k}_{T,2}",nBinsPt,minPt,maxPt,nBinsKt,0.,120.,nBinsKt,0.,120.);
 203   fOutput->Add(fh3PtTrigKt1Kt2FuCh);
 204
 205   fh3PtTrigDPhi1DPhi2Ch = new TH3F("fh3PtTrigDPhi1DPhi2Ch","fh3PtTrigDPhi1DPhi2Ch;#it{p}_{T,1} (GeV/#it{c});#it{k}_{T,1};#it{k}_{T,2}",nBinsPt,minPt,maxPt,36,-0.5*TMath::Pi(),1.5*TMath::Pi(),36,-0.5*TMath::Pi(),1.5*TMath::Pi());
 206   fOutput->Add(fh3PtTrigDPhi1DPhi2Ch);
 207
 208   fh3PtTrigDPhi1DPhi2FuCh = new TH3F("fh3PtTrigDPhi1DPhi2FuCh","fh3PtTrigDPhi1DPhi2FuCh;#it{p}_{T,1} (GeV/#it{c});#it{k}_{T,1};#it{k}_{T,2}",nBinsPt,minPt,maxPt,36,-0.5*TMath::Pi(),1.5*TMath::Pi(),36,-0.5*TMath::Pi(),1.5*TMath::Pi());
 209   fOutput->Add(fh3PtTrigDPhi1DPhi2FuCh);
 210
 211
 212   for(Int_t i=0; i<4; i++) {
 213     TString histoName = Form("fh3DiJetKtNEFPtAssoc_TrigBin%d",i);
 214     fh3DiJetKtNEFPtAssoc[i] = new TH3F(histoName.Data(),histoName.Data(),nBinsKt,0.,120.,50,0.,1.,nBinsPt,minPt,maxPt);
 215     fOutput->Add(fh3DiJetKtNEFPtAssoc[i]);
 216
 217     histoName = Form("fCentCorrPtAssocCh_TrigBin%d",i);
 218     fCentCorrPtAssocCh[i]  = new TH3F(histoName.Data(),histoName.Data(),10,0.,100.,10,0.,100.,nBinsPt,minPt,maxPt);
 219     fOutput->Add(fCentCorrPtAssocCh[i]);
 220
 221     histoName = Form("fCentCorrPtAssocFuCh_TrigBin%d",i);
 222     fCentCorrPtAssocFuCh[i]  = new TH3F(histoName.Data(),histoName.Data(),10,0.,100.,10,0.,100.,nBinsPt,minPt,maxPt);
 223     fOutput->Add(fCentCorrPtAssocFuCh[i]);
 224
 225     histoName = Form("fAjPtAssocCentCh_TrigBin%d",i);
 226     fAjPtAssocCentCh[i]  = new TH3F(histoName.Data(),histoName.Data(),50,0.,1.,nBinsPt,minPt,maxPt,20,0.,100.);
 227     fOutput->Add(fAjPtAssocCentCh[i]);
 228
 229     histoName = Form("fAjPtAssocCentFuCh_TrigBin%d",i);
 230     fAjPtAssocCentFuCh[i]  = new TH3F(histoName.Data(),histoName.Data(),50,0.,1.,nBinsPt,minPt,maxPt,20,0.,100.);
 231     fOutput->Add(fAjPtAssocCentFuCh[i]);
 232
 233     histoName = Form("fh3PtAssoc1PtAssoc2DPhi23Ch_TrigBin%d",i);
 234     fh3PtAssoc1PtAssoc2DPhi23Ch[i] = new TH3F(histoName.Data(),histoName.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt,nBinsDPhi,-0.5*TMath::Pi(),1.5*TMath::Pi());
 235     fOutput->Add(fh3PtAssoc1PtAssoc2DPhi23Ch[i]);
 236
 237     histoName = Form("fh3PtAssoc1PtAssoc2DPhi23FuCh_TrigBin%d",i);
 238     fh3PtAssoc1PtAssoc2DPhi23FuCh[i] = new TH3F(histoName.Data(),histoName.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt,nBinsDPhi,-0.5*TMath::Pi(),1.5*TMath::Pi());
 239     fOutput->Add(fh3PtAssoc1PtAssoc2DPhi23FuCh[i]);
 240   }
 241
 242   const Int_t nBinsSparseMatch = 7;
 243   const Int_t nBinsDPhiMatch = 80;
 244   const Int_t nBinsDEtaMatch = 80;
 245   const Int_t nBinsDR = 20;
 246   const Int_t nBinsFraction = 21;
 247   const Int_t nBinsType = 3;
 248   const Int_t nBinsMatch[nBinsSparseMatch] = {nBinsPt,nBinsPt,nBinsDPhiMatch,nBinsDEtaMatch,nBinsDR,nBinsFraction,nBinsType};
 249   //pTfull, pTch, deltaPhi, deltaEta, deltaR, fraction, jet type (leading,subleading,other)
 250   const Double_t xminMatch[nBinsSparseMatch]  = { minPt, minPt, -0.5,-0.5, 0. ,0.  ,0};
 251   const Double_t xmaxMatch[nBinsSparseMatch]  = { maxPt, maxPt,  0.5, 0.5, 0.5,1.05,3};
 252   if(fDoMatchFullCharged) {
 253     fhnMatchingFullCharged = new THnSparseF("fhnMatchingFullCharged","fhnMatchingFullCharged;#it{p}_{T,full} (GeV/#it{c});#it{p}_{T,ch} (GeV/#it{c});#Delta#varphi;#Delta#eta;#Delta R;f_{ch};type",
 254                                           nBinsSparseMatch,nBinsMatch,xminMatch,xmaxMatch);
 255     fOutput->Add(fhnMatchingFullCharged);
 256   }
 257
 258   // =========== Switch on Sumw2 for all histos ===========
 259   for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
 260     TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
 261     if (h1){
 262       h1->Sumw2();
 263       continue;
 264     }
 265     THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
 266     if(hn) {
 267       hn->Sumw2();
 268       continue;
 269     }
 270   }
 271
 272   TH1::AddDirectory(oldStatus);
 273
 274   PostData(1, fOutput); // Post data for ALL output slots > 0 here.
 275 }
 276
 277 //________________________________________________________________________
 278 Bool_t AliAnalysisTaskEmcalDiJetAna::FillHistograms()
 279 {
 280   // Fill histograms.
 281
 282   fh2CentRhoCh->Fill(fCent,fRhoChVal);
 283   fh2CentRhoScaled->Fill(fCent,fRhoFullVal);
 284
 285   Int_t nJetsFull = GetNJets(fContainerFull);
 286   for(Int_t ij=0; ij<nJetsFull; ij++) {
 287     AliEmcalJet *jet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ij, fContainerFull));
 288     if(!jet) continue; //jet not selected
 289
 290     Double_t jetPt = GetJetPt(jet,0);
 291     fh3PtEtaPhiJetFull->Fill(jetPt,jet->Eta(),jet->Phi());
 292   }
 293
 294   Int_t nJetsCh = GetNJets(fContainerCharged);
 295   for(Int_t ij=0; ij<nJetsCh; ij++) {
 296     AliEmcalJet *jet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ij, fContainerCharged));
 297     if(!jet) continue; //jet not selected
 298
 299       Double_t jetPt = GetJetPt(jet,1);
 300       fh3PtEtaPhiJetCharged->Fill(jetPt,jet->Eta(),jet->Phi());
 301   }
 302
 303   return kTRUE;
 304 }
 305
 306 //________________________________________________________________________
 307 Bool_t AliAnalysisTaskEmcalDiJetAna::Run()
 308 {
 309   // Run analysis code here, if needed. It will be executed before FillHistograms().
 310
 311   //Check if event is selected (vertex & pile-up)
 312   if(!SelectEvent())
 313     return kFALSE;
 314
 315   if(fRhoType==0) {
 316     fRhoFullVal = 0.;
 317     fRhoChVal = 0.;
 318   }
 319   if(fRhoType==1) {
 320     fRhoFullVal = GetRhoVal(fContainerFull);
 321     fRhoChVal = GetRhoVal(fContainerCharged);
 322   }
 323
 324   if(fDoFullFull)
 325     CorrelateJets(0);
 326
 327   // MatchFullAndChargedJets();
 328   if(fDoChargedCharged)   CorrelateJets(1);
 329
 330   if(fDoMatchFullCharged) FillMatchFullChargedHistos(fContainerFull,fContainerCharged);
 331
 332   if(fDoFullCharged)      {
 333     SetChargedFractionIndex();
 334     CorrelateJets(2);
 335   }
 336
 337   return kTRUE;  // If return kFALSE FillHistogram() will NOT be executed.
 338 }
 339
 340 //________________________________________________________________________
 341 void AliAnalysisTaskEmcalDiJetAna::CorrelateJets(const Int_t type) {
 342   //
 343   // Correlate jets and fill histos
 344   //
 345
 346   if( fJetCorrelationType==kCorrelateAll )
 347     CorrelateAllJets(type);
 348   else if( fJetCorrelationType==kCorrelateTwo )
 349     CorrelateTwoJets(type);
 350   else if( fJetCorrelationType==kCorrelateLS )
 351     CorrelateLeadingSubleadingJets(type);
 352
 353   return;
 354
 355 }
 356
 357 //________________________________________________________________________
 358 void AliAnalysisTaskEmcalDiJetAna::CorrelateTwoJets(const Int_t type) {
 359   //
 360   // Correlate jets and fill histos
 361   //
 362
 363   Int_t typet = 0;
 364   Int_t typea = 0;
 365   if(type==0) { //full-full
 366     typet = fContainerFull;
 367     typea = fContainerFull;
 368   }
 369   else if(type==1) { //charged-charged
 370     typet = fContainerCharged;
 371     typea = fContainerCharged;
 372   }
 373   else if(type==2) { //full-charged
 374     typet = fContainerFull;
 375     typea = fContainerCharged;
 376   }
 377   else {
 378     AliWarning(Form("%s: type %d of dijet correlation not defined!",GetName(),type));
 379     return;
 380   }
 381
 382   Int_t nJetsTrig  = GetNJets(typet);
 383   for(Int_t ijt=0; ijt<nJetsTrig; ijt++) {
 384
 385     AliEmcalJet *jetTrig = NULL;
 386     if(type==0) {
 387       jetTrig = static_cast<AliEmcalJet*>(GetJetFromArray(ijt, typet));
 388       if(TMath::Abs(jetTrig->Eta())>0.5)
 389         jetTrig = NULL;
 390     }
 391     else
 392       jetTrig = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ijt, typet));
 393
 394     if(!jetTrig)
 395       continue; //jet not selected
 396
 397     Double_t jetTrigPt = GetJetPt(jetTrig,typet);
 398
 399     if(jetTrigPt<fPtMinTriggerJet)
 400       continue;
 401
 402     AliEmcalJet *jetAssoc = GetLeadingAssociatedJet(typea,jetTrig);
 403     if(!jetAssoc)
 404       continue;
 405
 406     if(fDoPtBias) {
 407       if(type==0 || type==1) {
 408         if(GetJetPt(jetAssoc,typea)>jetTrigPt)
 409           continue;
 410       }
 411     }
 412
 413     FillDiJetHistos(jetTrig,jetAssoc, type);
 414
 415     //Look for second jet on away side - 3-jet events
 416     AliEmcalJet *jetAssoc2 = GetSecondLeadingAssociatedJet(typea,jetTrig);
 417     if(jetAssoc2)
 418       FillThreeJetHistos(jetTrig,jetAssoc,jetAssoc2,type);
 419
 420   }
 421 }
 422
 423 //________________________________________________________________________
 424 AliEmcalJet* AliAnalysisTaskEmcalDiJetAna::GetLeadingJet(const Int_t type) {
 425
 426   //Get associated jet which is the leading jet in the opposite hemisphere
 427
 428   Int_t cont = 0;
 429   if(type==0)  //full-full
 430     cont = fContainerFull;
 431   else if(type==1)  //charged-charged
 432     cont = fContainerCharged;
 433   else if(type==2)  //full-charged
 434     cont = fContainerFull;
 435
 436   Int_t nJets = GetNJets(cont);
 437   Double_t ptLead = -999;
 438   Int_t    iJetLead = -1;
 439   for(Int_t ij=0; ij<nJets; ij++) {
 440     AliEmcalJet *jet = NULL;
 441     if(type==0) {
 442       jet = static_cast<AliEmcalJet*>(GetJetFromArray(ij, cont));
 443       if(TMath::Abs(jet->Eta())>0.5)
 444         jet = NULL;
 445     }
 446     else
 447       jet = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ij, cont));
 448
 449     if(!jet)
 450       continue;
 451
 452     Double_t jetPt = GetJetPt(jet,cont);
 453
 454     if(jetPt>ptLead) {
 455       ptLead = jetPt;
 456       iJetLead = ij;
 457     }
 458   }
 459
 460   AliEmcalJet *jetLead = static_cast<AliEmcalJet*>(GetJetFromArray(iJetLead, cont));
 461
 462   return jetLead;
 463 }
 464
 465 //________________________________________________________________________
 466 AliEmcalJet* AliAnalysisTaskEmcalDiJetAna::GetLeadingAssociatedJet(const Int_t type, AliEmcalJet *jetTrig) {
 467
 468   //Get associated jet which is the leading jet in the opposite hemisphere
 469
 470   Int_t typea = 0;
 471   if(type==0)  //full-full
 472     typea = fContainerFull;
 473   else if(type==1)  //charged-charged
 474     typea = fContainerCharged;
 475   else if(type==2)  //full-charged
 476     typea = fContainerCharged;
 477
 478   AliEmcalJet *jetAssocLead = GetLeadingJetOppositeHemisphere(type, typea, jetTrig);
 479
 480   return jetAssocLead;
 481 }
 482
 483 //________________________________________________________________________
 484 AliEmcalJet* AliAnalysisTaskEmcalDiJetAna::GetSecondLeadingAssociatedJet(const Int_t type, AliEmcalJet *jetTrig) {
 485
 486   //Get associated jet which is the leading jet in the opposite hemisphere
 487
 488   Int_t typea = 0;
 489   if(type==0)  //full-full
 490     typea = fContainerFull;
 491   else if(type==1)  //charged-charged
 492     typea = fContainerCharged;
 493   else if(type==2)  //full-charged
 494     typea = fContainerCharged;
 495
 496   AliEmcalJet *jetAssocLead2 = GetSecondLeadingJetOppositeHemisphere(type, typea, jetTrig);
 497
 498   return jetAssocLead2;
 499 }
 500
 501 //________________________________________________________________________
 502 void AliAnalysisTaskEmcalDiJetAna::CorrelateAllJets(const Int_t type) {
 503   //
 504   // Correlate jets and fill histos
 505   //
 506
 507   Int_t typet = 0;
 508   Int_t typea = 0;
 509   if(type==0) { //full-full
 510     typet = fContainerFull;
 511     typea = fContainerFull;
 512   }
 513   else if(type==1) { //charged-charged
 514     typet = fContainerCharged;
 515     typea = fContainerCharged;
 516   }
 517   else if(type==2) { //full-charged
 518     typet = fContainerFull;
 519     typea = fContainerCharged;
 520   }
 521   else {
 522     AliWarning(Form("%s: type %d of dijet correlation not defined!",GetName(),type));
 523     return;
 524   }
 525
 526   Int_t nJetsTrig  = GetNJets(typet);
 527   Int_t nJetsAssoc = GetNJets(typea);
 528
 529   for(Int_t ijt=0; ijt<nJetsTrig; ijt++) {
 530     AliEmcalJet *jetTrig = NULL;
 531     if(type==0) {
 532       jetTrig = static_cast<AliEmcalJet*>(GetJetFromArray(ijt, typet));
 533       if(TMath::Abs(jetTrig->Eta())>0.5)
 534         jetTrig = NULL;
 535     }
 536     else
 537       jetTrig = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ijt, typet));
 538
 539     if(!jetTrig)
 540       continue; //jet not selected
 541
 542     Double_t jetTrigPt = GetJetPt(jetTrig,typet);
 543
 544     if(jetTrigPt<fPtMinTriggerJet)
 545       continue;
 546
 547     for(Int_t ija=0; ija<nJetsAssoc; ija++) {
 548       if(IsSameJet(ijt,ija,type)) continue;
 549
 550       AliEmcalJet *jetAssoc = NULL;
 551       if(type==0) {
 552         jetAssoc = static_cast<AliEmcalJet*>(GetJetFromArray(ija, typea));
 553         if(TMath::Abs(jetAssoc->Eta())>0.5)
 554           jetAssoc = NULL;
 555       }
 556       else
 557         jetAssoc = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ija, typea));
 558
 559       if(!jetAssoc)
 560         continue;
 561
 562       Double_t jetAssocPt = GetJetPt(jetAssoc,typea);
 563
 564       if(jetTrigPt>jetAssocPt)
 565         FillDiJetHistos(jetTrig,jetAssoc, type);
 566
 567     } // associate jet loop
 568   }//trigger jet loop
 569
 570 }
 571
 572 //________________________________________________________________________
 573 void AliAnalysisTaskEmcalDiJetAna::CorrelateLeadingSubleadingJets(const Int_t type) {
 574   //
 575   // Correlate leading jet in event with leading jet in opposite hemisphere
 576   //
 577
 578   Int_t typet = 0;
 579   Int_t typea = 0;
 580   if(type==0) { //full-full
 581     typet = fContainerFull;
 582     typea = fContainerFull;
 583   }
 584   else if(type==1) { //charged-charged
 585     typet = fContainerCharged;
 586     typea = fContainerCharged;
 587   }
 588   else if(type==2) { //full-charged
 589     typet = fContainerFull;
 590     typea = fContainerCharged;
 591   }
 592   else {
 593     AliWarning(Form("%s: type %d of dijet correlation not defined!",GetName(),type));
 594     return;
 595   }
 596
 597   AliEmcalJet *jetTrig = GetLeadingJet(typet);
 598   if(!jetTrig)
 599     return;
 600
 601   Double_t jetTrigPt = GetJetPt(jetTrig,typet);
 602
 603   if(jetTrigPt<fPtMinTriggerJet)
 604     return;
 605
 606   AliEmcalJet *jetAssoc = GetLeadingAssociatedJet(typea,jetTrig);
 607   if(!jetAssoc)
 608     return;
 609
 610   FillDiJetHistos(jetTrig,jetAssoc, type);
 611 }
 612
 613 //________________________________________________________________________
 614 void AliAnalysisTaskEmcalDiJetAna::FillDiJetHistos(const AliEmcalJet *jet1, const AliEmcalJet *jet2, const Int_t mode) {
 615   //
 616   // Fill histos
 617   // mode: full vs full        = 0
 618   //       charged vs charged  = 1
 619   //       full vs charged     = 2
 620   //
 621
 622   Int_t typet = 0;
 623   Int_t typea = 0;
 624   if(mode==0) { //full-full
 625     typet = fContainerFull;
 626     typea = fContainerFull;
 627   }
 628   else if(mode==1) { //charged-charged
 629     typet = fContainerCharged;
 630     typea = fContainerCharged;
 631   }
 632   else if(mode==2) { //full-charged
 633     typet = fContainerFull;
 634     typea = fContainerCharged;
 635   }
 636   else {
 637     AliWarning(Form("%s: mode %d of dijet correlation not defined!",GetName(),mode));
 638     return;
 639   }
 640
 641   Double_t jetTrigPt = GetJetPt(jet1,typet);
 642   Double_t jetAssocPt = GetJetPt(jet2,typea);
 643
 644   Double_t deltaPhi = GetDeltaPhi(jet1->Phi(),jet2->Phi());
 645
 646   Double_t kT = TMath::Abs(jetTrigPt*TMath::Sin(deltaPhi));
 647
 648   Double_t dijetEta = (jet1->Eta()+jet2->Eta())/2.;
 649
 650   Double_t aj = 0.;
 651   if((jetTrigPt+jetAssocPt)>0.) aj = (jetTrigPt-jetAssocPt)/(jetTrigPt+jetAssocPt);
 652
 653   Double_t diJetVars[7] = {jetTrigPt,jetAssocPt,deltaPhi,kT,dijetEta,fCent,aj};
 654
 655   if(mode==0)
 656     fhnDiJetVarsFull->Fill(diJetVars);
 657   else if(mode==1)
 658     fhnDiJetVarsCh->Fill(diJetVars);
 659   else if(mode==2)
 660     fhnDiJetVarsFullCharged->Fill(diJetVars);
 661
 662   Double_t dPhiMin = TMath::Pi() - 1./3.*TMath::Pi();
 663   Double_t dPhiMax = TMath::Pi() + 1./3.*TMath::Pi();
 664   Int_t trigBin = GetPtTriggerBin(jetTrigPt);
 665   if(mode==2) {
 666     if(trigBin>-1 && trigBin<4) {
 667       if(deltaPhi>dPhiMin && deltaPhi<dPhiMax)
 668         fh3DiJetKtNEFPtAssoc[trigBin]->Fill(kT, jet1->NEF(), jetAssocPt);
 669     }
 670   }
 671
 672   //Fill centrality correlation histos in case a dijet is present in acceptance
 673   Double_t centZNA = -1.;
 674   AliCentrality *aliCent = InputEvent()->GetCentrality();
 675   if (aliCent) {
 676     centZNA = aliCent->GetCentralityPercentile("ZNA");
 677     if(trigBin>-1 && trigBin<4) {
 678       if(deltaPhi>dPhiMin && deltaPhi<dPhiMax) {
 679         if(mode==1) {
 680           fCentCorrPtAssocCh[trigBin]->Fill(fCent,centZNA,jetAssocPt);
 681           fAjPtAssocCentCh[trigBin]->Fill(aj,jetAssocPt,fCent);
 682         }
 683         else if(mode==2) {
 684           fCentCorrPtAssocFuCh[trigBin]->Fill(fCent,centZNA,jetAssocPt);
 685           fAjPtAssocCentFuCh[trigBin]->Fill(aj,jetAssocPt,fCent);
 686         }
 687       }
 688     }
 689   }
 690
 691 }
 692
 693 //________________________________________________________________________
 694 void AliAnalysisTaskEmcalDiJetAna::FillThreeJetHistos(const AliEmcalJet *jet1, const AliEmcalJet *jet2, const AliEmcalJet *jet3, const Int_t mode) {
 695   //
 696   // Fill histos
 697   // mode: full vs full        = 0
 698   //       charged vs charged  = 1
 699   //       full vs charged     = 2
 700   //
 701
 702   Int_t typet = 0;
 703   Int_t typea = 0;
 704   if(mode==0) { //full-full
 705     typet = fContainerFull;
 706     typea = fContainerFull;
 707   }
 708   else if(mode==1) { //charged-charged
 709     typet = fContainerCharged;
 710     typea = fContainerCharged;
 711   }
 712   else if(mode==2) { //full-charged
 713     typet = fContainerFull;
 714     typea = fContainerCharged;
 715   }
 716   else {
 717     AliWarning(Form("%s: mode %d of dijet correlation not defined!",GetName(),mode));
 718     return;
 719   }
 720
 721   Double_t jetTrigPt = GetJetPt(jet1,typet);
 722   Double_t jetAssoc2Pt = GetJetPt(jet2,typea);
 723   Double_t jetAssoc3Pt = GetJetPt(jet3,typea);
 724
 725   Double_t deltaPhi12 = GetDeltaPhi(jet1->Phi(),jet2->Phi());
 726   Double_t deltaPhi13 = GetDeltaPhi(jet1->Phi(),jet3->Phi());
 727   Double_t deltaPhi23 = GetDeltaPhi(jet2->Phi(),jet3->Phi());
 728
 729   Double_t kT12 = TMath::Abs(jetTrigPt*TMath::Sin(deltaPhi12));
 730   Double_t kT13 = TMath::Abs(jetTrigPt*TMath::Sin(deltaPhi13));
 731
 732   Double_t dPhiMin = TMath::Pi() - 1./3.*TMath::Pi();
 733   Double_t dPhiMax = TMath::Pi() + 1./3.*TMath::Pi();
 734
 735   Int_t trigBin = GetPtTriggerBin(jetTrigPt);
 736
 737   if(jetAssoc2Pt>20. && jetAssoc3Pt>20.) {
 738     if(mode==1) {
 739       fh3PtTrigDPhi1DPhi2Ch->Fill(jetTrigPt,deltaPhi12,deltaPhi13);
 740       fh3PtAssoc1PtAssoc2DPhi23Ch[trigBin]->Fill(jetAssoc2Pt,jetAssoc3Pt,deltaPhi23);
 741     }
 742     else if(mode==1) {
 743       fh3PtTrigDPhi1DPhi2FuCh->Fill(jetTrigPt,deltaPhi12,deltaPhi13);
 744       fh3PtAssoc1PtAssoc2DPhi23FuCh[trigBin]->Fill(jetAssoc2Pt,jetAssoc3Pt,deltaPhi23);
 745     }
 746
 747     if(deltaPhi12>dPhiMin && deltaPhi12<dPhiMax) {
 748       if(mode==1)
 749         fh3PtTrigKt1Kt2Ch->Fill(jetTrigPt,kT12,kT13);
 750       else if(mode==2)
 751         fh3PtTrigKt1Kt2FuCh->Fill(jetTrigPt,kT12,kT13);
 752     }
 753   }
 754
 755 }
 756
 757 //________________________________________________________________________
 758 Int_t AliAnalysisTaskEmcalDiJetAna::GetPtTriggerBin(Double_t pt) {
 759
 760   Int_t binTrig = -1;
 761   if(pt>=20 && pt<40)
 762     binTrig = 0;
 763   else if(pt>=40 && pt<60)
 764     binTrig = 1;
 765   else if(pt>=60 && pt<80)
 766     binTrig = 2;
 767   else if(pt>=80 && pt<100)
 768     binTrig = 3;
 769
 770   return binTrig;
 771 }
 772
 773 //________________________________________________________________________
 774 void AliAnalysisTaskEmcalDiJetAna::FillMatchFullChargedHistos(Int_t cFull,Int_t cCharged) {
 775   //
 776   // Match full to charged jets and fill histo
 777   //
 778
 779   Int_t match = MatchFullAndChargedJets(cFull,cCharged);
 780   if(match==0) {
 781     AliDebug(11,Form("%s: matching failed",GetName()));
 782     return;
 783   }
 784
 785   for(int ig = 0;ig < GetNJets(cFull);++ig){
 786     AliEmcalJet *jetFull = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(ig, cFull));
 787     if(!jetFull) continue;
 788
 789     AliEmcalJet *jetCh = jetFull->ClosestJet();
 790     if(!jetCh) continue;
 791
 792     Double_t shFraction = GetFractionSharedPt(jetFull,jetCh);
 793     Double_t matchVars[7] = {
 794       jetFull->Pt(),
 795       jetCh->Pt(),
 796       GetDeltaPhi(jetFull->Phi(),jetCh->Phi()),
 797       jetFull->Eta()-jetCh->Eta(),
 798       GetDeltaR(jetFull,jetCh),
 799       shFraction,TMath::Min((Float_t)ig+0.5,2.5)
 800     };
 801       fhnMatchingFullCharged->Fill(matchVars);
 802
 803   }//loop over full jets
 804
 805 }
 806
 807
 808 //________________________________________________________________________
 809 Int_t AliAnalysisTaskEmcalDiJetAna::MatchFullAndChargedJets(Int_t cFull, Int_t cCharged) {
 810   //
 811   // Match charged jets to full jets
 812   //
 813
 814   if(GetNJets(cFull)<1) {
 815     AliDebug(2,Form("%s: no full jets: %d", GetName(),GetNJets(cFull)));
 816     return 0;
 817   }
 818
 819   if(GetNJets(cCharged)<1) {
 820     AliDebug(2,Form("%s: no charged jets: %d", GetName(),GetNJets(cCharged)));
 821     return 0;
 822   }
 823
 824   TClonesArray *cJetsFull = GetJetArray(cFull);
 825   TClonesArray *cJetsCharged = GetJetArray(cCharged);
 826
 827   if(!cJetsFull) {
 828     AliDebug(2,Form("%s: no full jet array",GetName()));
 829     return 0;
 830   }
 831
 832   if(!cJetsCharged) {
 833     AliDebug(2,Form("%s: no charged jet array",GetName()));
 834     return 0;
 835   }
 836
 837   if(!fMatchingDone) {
 838       MatchJetsGeo(cFull, cCharged, 0);
 839       return 1;
 840   } else {
 841     AliDebug(11,Form("%s: Matching already done before",GetName()));
 842     return 1;
 843   }
 844
 845 }
 846
 847 //_______________________________________________________________________
 848 void AliAnalysisTaskEmcalDiJetAna::Terminate(Option_t *)
 849 {
 850   // Called once at the end of the analysis.
 851 }
 852
 853
 854