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