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