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