PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskHJetDphi.cxx

   1 #include <TCanvas.h>
   2 #include <TChain.h>
   3 #include <TFormula.h>
   4 #include <TH1.h>
   5 #include <TH2.h>
   6 #include <TH3.h>
   7 #include <TProfile2D.h>
   8 #include <THnSparse.h>
   9 #include <TROOT.h>
  10 #include <TTree.h>
  11 #include <TArrayI.h>
  12 #include <TClonesArray.h>
  13 #include <TRandom3.h>
  14 #include <TFile.h>
  15 #include <TF1.h>
  16 #include <TLorentzVector.h>
  17 #include "TKey.h"
  18 #include "TList.h"
  19 #include "TSystem.h"
  20 #include "AliFJWrapper.h"
  21 #include "AliAODHandler.h"
  22 #include "AliAODEvent.h"
  23 #include "AliAODInputHandler.h"
  24 #include "AliAnalysisManager.h"
  25 #include "AliAnalysisTask.h"
  26 #include "AliCentrality.h"
  27 #include "AliAnalysisTaskHJetDphi.h"
  28 #include "AliESDEvent.h"
  29 #include "AliESDInputHandler.h"
  30 #include "AliESDtrack.h"
  31 #include "AliESDtrackCuts.h"
  32 #include "AliVParticle.h"
  33 #include "AliVTrack.h"
  34 #include "AliInputEventHandler.h"
  35 #include "AliMCEvent.h"
  36 #include "AliStack.h"
  37 #include "AliGenEventHeader.h"
  38 #include "AliGenPythiaEventHeader.h"
  39 #include "AliLog.h"
  40 #include "AliRhoParameter.h"
  41 #include "AliAODMCParticle.h"
  42 #include "AliNamedArrayI.h"
  43 #include "AliNamedString.h"
  44 #include "AliPicoTrack.h"
  45 #include "AliAnalysisTaskFastEmbedding.h"
  46 #include "AliEmcalJet.h"
  47 #include "AliAODJet.h"
  48 #include "AliAODJetEventBackground.h"
  49 #include "AliAnalysisHelperJetTasks.h"
  50
  51 #include <iostream>
  52 using std::cout;
  53 using std::cerr;
  54 using std::endl;
  55
  56 ClassImp(AliAnalysisTaskHJetDphi)
  57
  58 const Double_t pi = TMath::Pi();
  59 const Double_t kSector = pi/9;
  60
  61 //________________________________________________________________________
  62 AliAnalysisTaskHJetDphi::AliAnalysisTaskHJetDphi() :
  63   AliAnalysisTaskSE(),
  64   fVerbosity(0), fIsEmbedding(kFALSE), fAnaType(0), fPeriod("lhc11h"), fCollisionSystem("PbPb"),
  65   fIsMC(kFALSE), fAnalyzeMCTruth(kFALSE), fMC(0),
  66   fEvent(0x0), fESD(0x0), fAODIn(0x0), fAODOut(0x0), fAODExtension(0x0),
  67   fOfflineTrgMask(AliVEvent::kAny), fTriggerType(-1), fCentrality(-1), fMaxVtxZ(10),
  68   fEsdTrkCut(0x0), fEsdHybCut(0x0), fFilterMask(0), fRequireITSRefit(kTRUE), fRequireSharedClsCut(kTRUE),
  69   fIsInit(kFALSE), fNonStdFile(""), fMcParticleArrName(""), fMcParticleArray(0x0),  fMcParticlelMap(0x0),
  70   fEmbTrkArrName(""), fEmbTrkArray(0x0), fTrackArrName(""), fTrackArray(0x0),
  71   fTriggerTrkIndex(-1), fTriggerTrkPt(-1), fSwitchOnAvoidTpcHole(kFALSE), fAvoidTpcHole(0), fCutTPCBoundary(kFALSE), fDistToTPCBoundary(0.),
  72   fMinTrkPt(0.15), fMaxTrkPt(1e4), fMinTrkEta(-0.9), fMaxTrkEta(0.9), fMinTrkPhi(0), fMaxTrkPhi(2*pi),
  73   fRadius(0.4), fJetArrName(""), fPLJetArrName(""), fDLJetArrName(""), fJetArray(0x0), fPLJetArray(0x0), fDLJetArray(0x0),
  74   fRhoName(""), fRho(0x0), fRhoValue(0), fEvtBkg(0x0), fPtHardBinName(0x0), fPtHardBin(-1),
  75   fRunTrkQA(kFALSE), fRunJetQA(kFALSE), fRunSingleInclHJet(kFALSE),  fTTtype(0), fTTMinPt(9), fTTMaxPt(10), fJetPtMin(10),
  76   fRunPLHJet(kFALSE), fRunDLHJet(kFALSE), fRunLeadTrkQA(kFALSE), fStudyKtEffects(kFALSE), fKtValue(0), fRandom(0),
  77   fRunBkgFlow(kFALSE),
  78   fOutputList(0x0), fhEventStat(0x0), fhNTrials(0x0), fhPtHardBins(0x0)
  79 {
  80   // Constructor
  81
  82   // Define input and output slots here
  83   // Input slot #0 works with a TChain
  84   //DefineInput(0, TChain::Class());
  85   //DefineOutput(1, TList::Class());
  86   // Output slot #0 id reserved by the base class for AOD
  87
  88   for(Int_t i=0; i<4; i++)
  89     {
  90       fhVtxZ[i]                  = 0x0;
  91       fhCentrality[i]            = 0x0;
  92       fhRhoVsCent[i]             = 0x0;
  93
  94       fhTrkPt[i]                 = 0x0;
  95       fhTrkQA[i]                 = 0x0;
  96       fhTrkPtRes[i]              = 0x0;
  97       fhTrkPhiRes[i]             = 0x0;
  98
  99       fhNumberOfTT[i]            = 0x0;
 100       for(Int_t j=0; j<3; j++)
 101         {
 102           fhJetPt[i][j]          = 0x0;
 103           fhJetArea[i][j]        = 0x0;
 104           fhJetQA[i][j]          = 0x0;
 105
 106           fhTTPt[i][j]           = 0x0;
 107           fHJetPhiCorr[i][j]     = 0x0;
 108         }
 109       fHJetPhiCorrUp[i]          = 0x0;
 110       fHJetPhiCorrDown[i]        = 0x0;
 111
 112       fhLeadTrkQA[i]             = 0x0;
 113       fhKtEffects[i]             = 0x0;
 114     }
 115
 116   fAODfilterBits[0] = -1;
 117   fAODfilterBits[1] = -1;
 118 }
 119
 120 //________________________________________________________________________
 121 AliAnalysisTaskHJetDphi::AliAnalysisTaskHJetDphi(const char *name) :
 122   AliAnalysisTaskSE(name),
 123   fVerbosity(0), fIsEmbedding(kFALSE), fAnaType(0), fPeriod("lhc11h"), fCollisionSystem("PbPb"),
 124   fIsMC(kFALSE), fAnalyzeMCTruth(kFALSE), fMC(0),
 125   fEvent(0x0), fESD(0x0), fAODIn(0x0), fAODOut(0x0), fAODExtension(0x0),
 126   fOfflineTrgMask(AliVEvent::kAny), fTriggerType(-1), fCentrality(-1), fMaxVtxZ(10),
 127   fEsdTrkCut(0x0), fEsdHybCut(0x0), fFilterMask(0), fRequireITSRefit(kTRUE), fRequireSharedClsCut(kTRUE),
 128   fIsInit(kFALSE), fNonStdFile(""), fMcParticleArrName(""), fMcParticleArray(0x0),  fMcParticlelMap(0x0),
 129   fEmbTrkArrName(""), fEmbTrkArray(0x0), fTrackArrName(""), fTrackArray(0x0),
 130   fTriggerTrkIndex(-1), fTriggerTrkPt(-1), fSwitchOnAvoidTpcHole(kFALSE), fAvoidTpcHole(0), fCutTPCBoundary(kFALSE), fDistToTPCBoundary(0.),
 131   fMinTrkPt(0.15), fMaxTrkPt(1e4), fMinTrkEta(-0.9), fMaxTrkEta(0.9), fMinTrkPhi(0), fMaxTrkPhi(2*pi),
 132   fRadius(0.4), fJetArrName(""), fPLJetArrName(""), fDLJetArrName(""), fJetArray(0x0), fPLJetArray(0x0), fDLJetArray(0x0),
 133   fRhoName(""), fRho(0x0), fRhoValue(0), fEvtBkg(0x0), fPtHardBinName(0x0), fPtHardBin(-1),
 134   fRunTrkQA(kFALSE), fRunJetQA(kFALSE), fRunSingleInclHJet(kFALSE),  fTTtype(0), fTTMinPt(9), fTTMaxPt(10), fJetPtMin(10),
 135   fRunPLHJet(kFALSE), fRunDLHJet(kFALSE), fRunLeadTrkQA(kFALSE), fStudyKtEffects(kFALSE), fKtValue(0), fRandom(0),
 136   fRunBkgFlow(kFALSE),
 137   fOutputList(0x0), fhEventStat(0x0), fhNTrials(0x0), fhPtHardBins(0x0)
 138 {
 139   // Constructor
 140
 141   // Define input and output slots here
 142   // Input slot #0 works with a TChain
 143   DefineInput(0, TChain::Class());
 144   DefineOutput(1, TList::Class());
 145   // Output slot #0 id reserved by the base class for AOD
 146
 147   for(Int_t i=0; i<4; i++)
 148     {
 149       fhVtxZ[i]                  = 0x0;
 150       fhCentrality[i]            = 0x0;
 151       fhRhoVsCent[i]             = 0x0;
 152
 153       fhTrkPt[i]                 = 0x0;
 154       fhTrkQA[i]                 = 0x0;
 155       fhTrkPtRes[i]              = 0x0;
 156       fhTrkPhiRes[i]             = 0x0;
 157
 158       fhNumberOfTT[i]            = 0x0;
 159       for(Int_t j=0; j<3; j++)
 160         {
 161           fhJetPt[i][j]          = 0x0;
 162           fhJetArea[i][j]        = 0x0;
 163           fhJetQA[i][j]          = 0x0;
 164
 165           fhTTPt[i][j]           = 0x0;
 166           fHJetPhiCorr[i][j]     = 0x0;
 167         }
 168
 169       fHJetPhiCorrUp[i]          = 0x0;
 170       fHJetPhiCorrDown[i]        = 0x0;
 171
 172       fhLeadTrkQA[i]             = 0x0;
 173       fhKtEffects[i]             = 0x0;
 174     }
 175
 176   fAODfilterBits[0] = -1;
 177   fAODfilterBits[1] = -1;
 178 }
 179
 180
 181 //________________________________________________________________________
 182 AliAnalysisTaskHJetDphi::~AliAnalysisTaskHJetDphi()
 183 {
 184   //Destructor
 185   if(fRandom)      delete fRandom;
 186   if(fEsdTrkCut)   delete fEsdTrkCut;
 187   if(fEsdHybCut)   delete fEsdHybCut;
 188   if(fOutputList) { fOutputList->Delete(); delete fOutputList;}
 189 }
 190
 191 //________________________________________________________________________
 192 void AliAnalysisTaskHJetDphi::UserCreateOutputObjects()
 193 {
 194   // Create histograms
 195
 196   const Int_t nTrkPtBins = 100;
 197   const Float_t lowTrkPtBin=0, upTrkPtBin=100;
 198   const Int_t nJetPtBins = 300;
 199   const Float_t lowJetPtBin=-100, upJetPtBin=200;
 200
 201   // track QA
 202   const Int_t dimTrkqa = 4;
 203   const Int_t nBinsTrkqa[dimTrkqa]     = {nTrkPtBins/5,  36,  40, 10};
 204   const Double_t lowBinTrkqa[dimTrkqa] = {lowTrkPtBin,   0,  -1,  0};
 205   const Double_t hiBinTrkqa[dimTrkqa]  = {upTrkPtBin,    360, 1,  10};
 206
 207   const Int_t dimTrkRes = 5;
 208   const Int_t nBinsTrkRes[dimTrkRes]     = {nTrkPtBins,    50,  50,  3, 10};
 209   const Double_t lowBinTrkRes[dimTrkRes] = {lowTrkPtBin,   0,   0,   0,  0};
 210   const Double_t hiBinTrkRes[dimTrkRes]  = {upTrkPtBin,    0.5, 0.5, 3,  10};
 211
 212   const Int_t dimPhiRes = 4;
 213   const Int_t nBinsPhiRes[dimPhiRes]     = {nTrkPtBins,    200,     3, 10};
 214   const Double_t lowBinPhiRes[dimPhiRes] = {lowTrkPtBin,   -0.00995, 0, 0};
 215   const Double_t hiBinPhiRes[dimPhiRes]  = {upTrkPtBin,    0.01005,  3, 10};
 216
 217   // jet QA
 218   const Int_t dimJetpt = 4;
 219   const Int_t nBinsJetpt[dimJetpt]     = {nJetPtBins,    300, 10, 10};
 220   const Double_t lowBinJetpt[dimJetpt] = {lowJetPtBin,   0,   0,  0};
 221   const Double_t hiBinJetpt[dimJetpt]  = {upJetPtBin,    300, 10, 10};
 222
 223   const Int_t dimJetA = 4;
 224   const Int_t nBinsJetA[dimJetA]     = {nJetPtBins,    100, 10, 10};
 225   const Double_t lowBinJetA[dimJetA] = {lowJetPtBin,   0,   0,  0};
 226   const Double_t hiBinJetA[dimJetA]  = {upJetPtBin,    1,   10, 10};
 227
 228   const Int_t dimJetqa = 7;
 229   const Int_t nBinsJetqa[dimJetqa]     = {nJetPtBins/5, 36,  24,  6,   100, 10, 11};
 230   const Double_t lowBinJetqa[dimJetqa] = {lowJetPtBin,   0,  -0.6, 0,   0,   0,  0};
 231   const Double_t hiBinJetqa[dimJetqa]  = {upJetPtBin,    360, 0.6, 1.2, 500, 10, 11};
 232
 233   // h-jet analysis
 234   const Int_t dimTT = 3;
 235   const Int_t nBinsTT[dimTT]     = {nTrkPtBins,  10,  11};
 236   const Double_t lowBinTT[dimTT] = {lowTrkPtBin, 0,   0};
 237   const Double_t hiBinTT[dimTT]  = {upTrkPtBin,  100, 11};
 238
 239   const Int_t dimCor = 8;
 240   const Int_t nBinsCor[dimCor]     = {nTrkPtBins, nJetPtBins,  140,     6,   10, 40,    11, 10};
 241   const Double_t lowBinCor[dimCor] = {lowTrkPtBin,lowJetPtBin, pi-4.95, 0,   0,  -1.95, 0,  0};
 242   const Double_t hiBinCor[dimCor]  = {upTrkPtBin, upJetPtBin,  pi+2.05, 1.2, 100, 2.05, 11, 10};
 243
 244   // Leading track QA
 245   const Int_t dimLeadTrkqa = 5;
 246   const Int_t nBinsLeadTrkqa[dimLeadTrkqa]     = {nTrkPtBins,  200,  80,   55, 10};
 247   const Double_t lowBinLeadTrkqa[dimLeadTrkqa] = {lowTrkPtBin,   0,  -0.4, 0,   0};
 248   const Double_t hiBinLeadTrkqa[dimLeadTrkqa]  = {upTrkPtBin,  200,  0.4,  1.1, 100};
 249
 250   // kt effects
 251   const Int_t dimKt = 5;
 252   const Int_t nBinsKt[dimKt]     = {nTrkPtBins,   20,  81,    10, 11};
 253   const Double_t lowBinKt[dimKt] = {lowTrkPtBin,   0,  -40.5, 0,  0};
 254   const Double_t hiBinKt[dimKt]  = {upTrkPtBin,  200,  40.5,  10, 11};
 255
 256    // Called once
 257    OpenFile(1);
 258    if(!fOutputList) fOutputList = new TList;
 259    fOutputList->SetOwner(kTRUE);
 260
 261    Bool_t oldStatus = TH1::AddDirectoryStatus();
 262    if(fAnaType==1) TH1::AddDirectory(kFALSE);
 263
 264   fhEventStat = new TH1F("fhEventStat","Event statistics for jet analysis",8,0,8);
 265   fhEventStat->GetXaxis()->SetBinLabel(1,"All");
 266   fhEventStat->GetXaxis()->SetBinLabel(2,"PS");
 267   fhEventStat->GetXaxis()->SetBinLabel(3,"Vtx");
 268   fhEventStat->GetXaxis()->SetBinLabel(4,"Vtx+10cm");
 269   fhEventStat->GetXaxis()->SetBinLabel(5,"kMB");
 270   fhEventStat->GetXaxis()->SetBinLabel(6,"kCentral");
 271   fhEventStat->GetXaxis()->SetBinLabel(7,"kSemiCentral");
 272   fhEventStat->GetXaxis()->SetBinLabel(8,"kJetService");
 273   fOutputList->Add(fhEventStat);
 274
 275   fhNTrials = new TH1F("fhNTrials","Number of trials",1,0,1);
 276   fOutputList->Add(fhNTrials);
 277
 278   fhPtHardBins = new TH1F("fhPtHardBins","Number of events in each pT hard bin",11,0,11);
 279   fOutputList->Add(fhPtHardBins);
 280
 281   const char *triggerName[4] = {"kMB","kEGA","kEJE","MC"};
 282   const char *dataType[3]    = {"", "_DL","_PL"};
 283   const char *dataName[3]    = {"Data","DL","PL"};
 284
 285   Double_t newbins[7] = {0,0.07,0.2,0.4,0.6,0.8,1};
 286
 287   for(Int_t i=0; i<4; i++)
 288     {
 289       if( fAnalyzeMCTruth )
 290         {
 291           if(i!=3) continue;
 292         }
 293       else
 294         {
 295           if( fPeriod=="lhc11a" && i>1 ) continue;
 296           if( fPeriod=="lhc10h" && i!=0 ) continue;
 297           if( fPeriod=="lhc11h" && i!=0 ) continue;
 298           if( fPeriod.Contains("lhc12a15a") && i!=0 ) continue;
 299           if( fPeriod.Contains("lhc12a15e") && i!=0 ) continue;
 300         }
 301
 302       fhVtxZ[i] = new TH1F(Form("%s_fhVtxZ",triggerName[i]),Form("%s: z distribution of event vertexz;z(cm)",triggerName[i]),400,-20,20);
 303       fOutputList->Add(fhVtxZ[i]);
 304
 305       fhCentrality[i] = new TH1F(Form("%s_fhCentrality",triggerName[i]),Form("%s: Event centrality;centrality",triggerName[i]),100,0,100);
 306       fOutputList->Add(fhCentrality[i]);
 307
 308       fhRhoVsCent[i] = new TH2F(Form("%s_fhRhoVsCent",triggerName[i]),Form("%s: Rho vs centrality (R=%1.1f);centrality;Rho",triggerName[i],fRadius),100,0,100,300,0,300);
 309       fOutputList->Add(fhRhoVsCent[i]);
 310
 311       if(fRunTrkQA)
 312         {
 313           fhTrkPt[i] = new TH2F(Form("%s_fhTrkPt",triggerName[i]),Form("%s: Track p_{T} vs centrality;p_{T}^{track} (GeV/c);Centrality",triggerName[i]),nTrkPtBins,lowTrkPtBin,upTrkPtBin,10,0,100);
 314           fOutputList->Add(fhTrkPt[i]);
 315
 316           fhTrkQA[i] = new THnSparseF(Form("%s_fhTrkQA",triggerName[i]),Form("%s: track p_{T} vs #phi vs #eta vs centrality;p_{T,track} (GeV/c);#phi;#eta;centrality",triggerName[i]),dimTrkqa,nBinsTrkqa,lowBinTrkqa,hiBinTrkqa);
 317           fOutputList->Add(fhTrkQA[i]);
 318
 319           fhTrkPtRes[i] = new THnSparseF(Form("%s_fhTrkPtRes",triggerName[i]),Form("%s: track p_{T} vs resolution vs (p_{T}^{gen}-p_{T}^{rec})/p_{T}^{gen} vs type vs centrality;p_{T,track} (GeV/c);#sigma(p_{T})/p_{T};type;centrality",triggerName[i]),dimTrkRes,nBinsTrkRes,lowBinTrkRes,hiBinTrkRes);
 320           fOutputList->Add(fhTrkPtRes[i]);
 321
 322           fhTrkPhiRes[i] = new THnSparseF(Form("%s_fhTrkPhiRes",triggerName[i]),Form("%s: track p_{T} vs #varphi^{gen}-#varphi^{rec} vs type vs centrality;p_{T,track} (GeV/c);#Delta#varphi;type;centrality",triggerName[i]),dimPhiRes,nBinsPhiRes,lowBinPhiRes,hiBinPhiRes);
 323           fOutputList->Add(fhTrkPhiRes[i]);
 324         }
 325
 326       for(Int_t j=0; j<3; j++)
 327         {
 328           if(!fRunDLHJet && j==1) continue;
 329           if(!fRunPLHJet && j==2) continue;
 330           if(fRunJetQA)
 331             {
 332               fhJetPt[i][j] = new THnSparseF(Form("%s_fhJetPt%s",triggerName[i],dataType[j]),Form("%s-%s: jet p_{T} vs raw jet p_{T} vs centrality vs pt hard bin (R=%1.1f);p_{T,jet}^{ch} (GeV/c);p_{T,jet}^{raw} (GeV/c);centrality;ptHardBin",dataName[j],triggerName[i],fRadius),dimJetpt,nBinsJetpt,lowBinJetpt,hiBinJetpt);
 333               fOutputList->Add(fhJetPt[i][j]);
 334
 335               fhJetArea[i][j] = new THnSparseF(Form("%s_fhJetArea%s",triggerName[i],dataType[j]),Form("%s-%s: jet p_{T} vs area vs centrality vs pt hard bin (R=%1.1f);p_{T,jet}^{ch} (GeV/c);area;centrality;ptHardBin",dataName[j],triggerName[i],fRadius),dimJetA,nBinsJetA,lowBinJetA,hiBinJetA);
 336               fOutputList->Add(fhJetArea[i][j]);
 337
 338               fhJetQA[i][j] = new THnSparseF(Form("%s_fhJetQA%s",triggerName[i],dataType[j]),Form("%s-%s: jet p_{T} vs #phi vs #eta vs area vs # of constituents vs centrality vs pt hard bin (R=%1.1f);p_{T,jet}^{ch} (GeV/c);#phi;#eta;area;# of constituents;centrality;ptHardBin",dataName[j],triggerName[i],fRadius),dimJetqa,nBinsJetqa,lowBinJetqa,hiBinJetqa);
 339               fhJetQA[i][j]->SetBinEdges(3,newbins);
 340               fOutputList->Add(fhJetQA[i][j]);
 341             }
 342
 343           if(fRunSingleInclHJet)
 344             {
 345               fhTTPt[i][j] = new THnSparseF(Form("%s_fhTTPt%s",triggerName[i],dataType[j]),Form("%s-%s: TT p_{T} vs centrality vs pT hard bin;p_{T,TT}^{ch} (GeV/c);centrality;ptHardBin",dataName[j],triggerName[i]),dimTT,nBinsTT,lowBinTT,hiBinTT);
 346               fOutputList->Add(fhTTPt[i][j]);
 347
 348               fHJetPhiCorr[i][j] = new THnSparseF(Form("%s_fHJetPhiCorr%s",triggerName[i],dataType[j]),Form("%s-%s: single inclusive TT p_{T} vs recoil jet p_{T} vs #Delta#varphi vs area vs centrality vs #Delta#eta vs pt hard bin vs event bin (R=%1.1f);TT p_{T} (GeV/c);p_{T,jet}^{ch} (GeV/c);#Delta#varphi;area;centrality;#Delta#eta;ptHardBin;EventBin",dataName[j],triggerName[i],fRadius),dimCor,nBinsCor,lowBinCor,hiBinCor);
 349               fHJetPhiCorr[i][j]->SetBinEdges(3,newbins);
 350               fOutputList->Add(fHJetPhiCorr[i][j]);
 351             }
 352         }
 353
 354
 355       if(fRunSingleInclHJet)
 356         {
 357           fhNumberOfTT[i] = new TH1F(Form("%s_fhNumberOfTT",triggerName[i]), Form("%s: number of TT",triggerName[i]),6,0,6);
 358           fOutputList->Add(fhNumberOfTT[i]);
 359
 360
 361           if(fRunBkgFlow)
 362             {
 363               fHJetPhiCorrUp[i] = new THnSparseF(Form("%s_fHJetPhiCorrUp",triggerName[i]),Form("%s: single inclusive TT p_{T} vs recoil jet p_{T} vs #Delta#varphi vs area vs centrality vs #Delta#eta vs pt hard bin vs event bin (R=%1.1f);TT p_{T} (GeV/c);p_{T,jet}^{ch} (GeV/c);#Delta#varphi;area;centrality;#Delta#eta;ptHardBin;EventBin",triggerName[i],fRadius),dimCor,nBinsCor,lowBinCor,hiBinCor);
 364               fOutputList->Add(fHJetPhiCorrUp[i]);
 365
 366               fHJetPhiCorrDown[i] = new THnSparseF(Form("%s_fHJetPhiCorrDown",triggerName[i]),Form("%s: single inclusive TT p_{T} vs recoil jet p_{T} vs #Delta#varphi vs area vs centrality vs #Delta#eta vs pt hard bin vs event bin (R=%1.1f);TT p_{T} (GeV/c);p_{T,jet}^{ch} (GeV/c);#Delta#varphi;area;centrality;#Delta#eta;ptHardBin;EventBin",triggerName[i],fRadius),dimCor,nBinsCor,lowBinCor,hiBinCor);
 367               fOutputList->Add(fHJetPhiCorrDown[i]);
 368             }
 369         }
 370
 371       if(fRunLeadTrkQA)
 372         {
 373           fhLeadTrkQA[i] = new THnSparseF(Form("%s_fhLeadTrkQA",triggerName[i]),Form("%s: p_{T,trk}^{leading} vs p_{T,jet}^{full} vs #Delta#varphi vs z vs centrality;p_{T,trk}^{leading} (GeV/c); p_{T,jet}^{full} (GeV/c);#Delta#varphi;z;centrality",triggerName[i]),dimLeadTrkqa,nBinsLeadTrkqa,lowBinLeadTrkqa,hiBinLeadTrkqa);
 374           fOutputList->Add(fhLeadTrkQA[i]);
 375         }
 376
 377       if(fStudyKtEffects)
 378         {
 379           fhKtEffects[i] = new THnSparseF(Form("%s_fhKtEffects",triggerName[i]),Form("%s: TT p_{T} vs recoil jet p_{T} vs k_{t} vs centrality vs pt hard bin (R=%1.1f);TT p_{T} (GeV/c);p_{T,jet}^{ch} (GeV/c);k_{t} (GeV/c);centrality;ptHardBin",triggerName[i],fRadius),dimKt,nBinsKt,lowBinKt,hiBinKt);
 380           fOutputList->Add(fhKtEffects[i]);
 381         }
 382     }
 383
 384   //error calculation in THnSparse
 385   Int_t nObj = fOutputList->GetEntries();
 386   for(Int_t i=0; i<nObj; i++)
 387     {
 388       TObject *obj = (TObject*) fOutputList->At(i);
 389       if (obj->IsA()->InheritsFrom( "THnSparse" ))
 390         {
 391           THnSparseF *hn = (THnSparseF*)obj;
 392           hn->Sumw2();
 393         }
 394     }
 395
 396   if(fRunTrkQA)
 397     {
 398       fEsdTrkCut = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kFALSE,1);
 399       fEsdTrkCut->SetMaxDCAToVertexXY(2.4);
 400       fEsdTrkCut->SetMaxDCAToVertexZ(3.2);
 401       fEsdTrkCut->SetDCAToVertex2D(kTRUE);
 402       fEsdTrkCut->SetMaxChi2TPCConstrainedGlobal(36);
 403       fEsdTrkCut->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
 404
 405       fEsdHybCut = new AliESDtrackCuts(*fEsdTrkCut);
 406       fEsdHybCut->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kNone);
 407     }
 408
 409   fRandom = new TRandom3(0);
 410
 411   PrintConfig();
 412
 413   if(fAnaType==1) TH1::AddDirectory(oldStatus);
 414   PostData(1, fOutputList);
 415 }
 416
 417 //________________________________________________________________________
 418 Bool_t AliAnalysisTaskHJetDphi::UserNotify()
 419 {
 420   AliInfo("User Nofity");
 421
 422   Int_t runNumber = InputEvent()->GetRunNumber();
 423
 424   fAvoidTpcHole = 0;
 425   if(fSwitchOnAvoidTpcHole)
 426     {
 427       Int_t runs_iroc[28] = {169975, 169981, 170038, 170040, 170083, 170084, 170085, 170088, 170089, 170091, 170152, 170155, 170159, 170163, 170193, 170195, 170203, 170204, 170205, 170228, 170230, 170264, 170268, 170269, 170270, 170306, 170308, 170309};
 428       Int_t runs_oroc[23] = {169591, 169590, 169588, 169587, 169586, 169584, 169557, 169555, 169554, 169553, 169550, 169515, 169512, 169506, 169504, 169498, 169475, 169420, 169418, 169099, 169040, 169045, 169044};
 429
 430       for(Int_t i=0; i<28; i++)
 431         {
 432           if(runNumber==runs_iroc[i])
 433             {
 434               fAvoidTpcHole = 1;
 435               break;
 436             }
 437         }
 438       for(Int_t i=0; i<23; i++)
 439         {
 440           if(runNumber==runs_oroc[i])
 441             {
 442               fAvoidTpcHole = 2;
 443               break;
 444             }
 445         }
 446     }
 447
 448   if(fIsMC)
 449     {
 450       TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 451       TFile *currFile = tree->GetCurrentFile();
 452       TString fileName(currFile->GetName());
 453       if(fileName.Contains("root_archive.zip#"))
 454         {
 455           Ssiz_t pos = fileName.Index("#",0,TString::kExact);
 456           fileName.Replace(pos+1,20,"");
 457         }
 458       else
 459         {
 460           fileName.ReplaceAll(gSystem->BaseName(fileName.Data()),"");
 461         }
 462
 463       TFile *fxsec = TFile::Open(Form("%s%s",fileName.Data(),"pyxsec_hists.root"),"read");
 464       if(fxsec)
 465         {
 466           TKey *key = (TKey*)fxsec->GetListOfKeys()->At(0);
 467           if(key)
 468             {
 469               TList *list = dynamic_cast<TList*>(key->ReadObj());
 470               if(list)
 471                 {
 472                   fhNTrials->Fill(0.5, ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1));
 473                 }
 474               else
 475                 return kFALSE;
 476             }
 477           else
 478             return kFALSE;
 479         }
 480       else
 481         {
 482           fxsec = TFile::Open(Form("%s%s",fileName.Data(),"pyxsec.root"),"read");
 483           TTree *xtree = (TTree*)fxsec->Get("Xsection");
 484           if(xtree)
 485             {
 486               UInt_t ntrials = 0;
 487               xtree->SetBranchAddress("ntrials",&ntrials);
 488               xtree->GetEntry(0);
 489               fhNTrials->Fill(0.5, ntrials);
 490             }
 491           else
 492             return kFALSE;
 493         }
 494       fxsec->Close();
 495     }
 496   return kTRUE;
 497 }
 498
 499 //________________________________________________________________________
 500 void AliAnalysisTaskHJetDphi::UserExec(Option_t *)
 501 {
 502   // Main loop, called for each event.
 503
 504   fTriggerType = -1;
 505
 506   // Get pointers to input events
 507   fEvent = InputEvent();
 508   if (!fEvent)
 509     {
 510       AliError("Input event not available");
 511       return;
 512     }
 513
 514   if(fIsMC)
 515     {
 516       fMC  = MCEvent();
 517       if (!fMC)
 518         {
 519           AliError("MC event available");
 520           return;
 521         }
 522     }
 523
 524   fESD=dynamic_cast<AliESDEvent*>(InputEvent());
 525   if (!fESD)
 526     {
 527       fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
 528     }
 529   else
 530     {
 531       if(fAnaType==1) fEvent = AODEvent();
 532     }
 533   if(fAnaType==1)
 534     {
 535       fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
 536       if(fNonStdFile.Length()!=0)
 537         {
 538           // case that we have an AOD extension we need can fetch the jets from the extended output
 539           AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 540           fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
 541           if(!fAODExtension)
 542             {
 543               if(fVerbosity>1) printf("AODExtension not found for %s",fNonStdFile.Data());
 544             }
 545         }
 546     }
 547
 548   fhEventStat->Fill(0.5);
 549
 550   if(fVerbosity>1)
 551     {
 552       TList *list = 0x0;
 553       if(fAnaType==0) list = fEvent->GetList();
 554       else            list = fAODOut->GetList();
 555       for(Int_t i=0; i<list->GetEntries(); i++)
 556         {
 557           TObject *obj = (TObject*)list->At(i);
 558           cout<<i<<": "<<obj->GetName()<<" : "<<obj->ClassName()<<endl;
 559         }
 560     }
 561   // Retrieve arraies from memory
 562   if(!fIsInit) fIsInit = RetrieveArraies();
 563   if(!fIsInit) return;
 564
 565   if(fIsEmbedding)
 566     {
 567       if(fAnaType==0)
 568         {
 569           TString fileName = fPtHardBinName->GetString();
 570           fileName.Remove(0,51);
 571           fileName.Remove(fileName.Index("/"));
 572           fPtHardBin = fileName.Atoi();
 573         }
 574       if(fAnaType==1)
 575         {
 576           Double_t pthard = AliAnalysisTaskFastEmbedding::GetPtHard();
 577           fPtHardBin = GetPtHardBin(pthard);
 578         }
 579     }
 580
 581   // physics selection
 582   UInt_t trigger = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 583   if(fOfflineTrgMask==AliVEvent::kAny) fTriggerType = 0;
 584   else
 585     {
 586       if(trigger & fOfflineTrgMask) fTriggerType=0;
 587
 588       if(fPeriod.Contains("lhc11a",TString::kIgnoreCase))
 589         {
 590           if (trigger & AliVEvent::kMB)     { fTriggerType=0; }
 591           if (trigger & AliVEvent::kEMC1)   { fTriggerType=1; }
 592         }
 593     }
 594   if(fTriggerType==-1) return;
 595   if(fAnalyzeMCTruth) fTriggerType = 3;
 596   fhEventStat->Fill(1.5);
 597
 598   // Vertex cut
 599   const AliVVertex* vtx = fEvent->GetPrimaryVertex();
 600   if (!vtx || vtx->GetNContributors()<1) return;
 601   fhEventStat->Fill(2.5);
 602   fhVtxZ[fTriggerType]->Fill(vtx->GetZ());
 603   if (TMath::Abs(vtx->GetZ())>fMaxVtxZ)  return;
 604   fhEventStat->Fill(3.5);
 605
 606   if(fTriggerType==0)
 607     {
 608       if(trigger & AliVEvent::kCentral) fhEventStat->Fill(5.5);
 609       else if (trigger & AliVEvent::kCentral) fhEventStat->Fill(6.5);
 610       else fhEventStat->Fill(4.5);
 611     }
 612
 613   if(!AliAnalysisHelperJetTasks::Selected()) return;
 614   fhEventStat->Fill(7.5);
 615
 616     // GetCentrality
 617   if(fCollisionSystem=="PbPb")
 618     {
 619       AliCentrality *centrality = fEvent->GetCentrality();
 620       if (centrality) fCentrality = centrality->GetCentralityPercentile("V0M");
 621       else            fCentrality = 99;
 622     }
 623   else if(fCollisionSystem=="pp")
 624     fCentrality = 0;
 625   fhCentrality[fTriggerType]->Fill(fCentrality);
 626
 627   // Get Rho value
 628   if(fCollisionSystem=="PbPb")
 629     {
 630       if(fAnaType==0) fRhoValue = fRho->GetVal();
 631       if(fAnaType==1) fRhoValue = fEvtBkg->GetBackground(0);
 632     }
 633   else if(fCollisionSystem=="pp")
 634     {
 635       fRhoValue = 0;
 636     }
 637   fhRhoVsCent[fTriggerType]->Fill(fCentrality,fRhoValue);
 638   fhPtHardBins->Fill(fPtHardBin);
 639
 640   if(fRunSingleInclHJet)
 641     {
 642       Double_t trigPt = -1, trigPhi = -999, trigEta = -999;
 643       Int_t nTrig = FindSingleIncTrigger(fTrackArray, trigPt, trigPhi, trigEta, fIsEmbedding);
 644       if(nTrig>0) fhNumberOfTT[fTriggerType]->Fill(nTrig);
 645       RunSingleInclHJetCorr(trigPt, trigPhi, trigEta, fJetArray, fRhoValue, fhTTPt[fTriggerType][0], fHJetPhiCorr[fTriggerType][0]);
 646       if(fRunBkgFlow)
 647         {
 648           RunSingleInclHJetCorr(trigPt, trigPhi, trigEta, fJetArray, fRhoValue+1, 0x0, fHJetPhiCorrUp[fTriggerType]);
 649           RunSingleInclHJetCorr(trigPt, trigPhi, trigEta, fJetArray, fRhoValue-1, 0x0, fHJetPhiCorrDown[fTriggerType]);
 650         }
 651
 652       if(fIsEmbedding)
 653         {
 654           if(fRunDLHJet)
 655             {
 656               FindSingleIncTrigger(fTrackArray, trigPt, trigPhi, trigEta, 1);
 657               RunSingleInclHJetCorr(trigPt, trigPhi, trigEta, fDLJetArray, 0, fhTTPt[fTriggerType][1], fHJetPhiCorr[fTriggerType][1]);
 658             }
 659
 660           if(fRunPLHJet)
 661             {
 662               FindSingleIncTrigger(fMcParticleArray, trigPt, trigPhi, trigEta, 2);
 663               RunSingleInclHJetCorr(trigPt, trigPhi, trigEta, fPLJetArray, 0, fhTTPt[fTriggerType][2], fHJetPhiCorr[fTriggerType][2]);
 664             }
 665         }
 666     }
 667
 668   if(fRunJetQA)
 669     {
 670       RunJetQA(fJetArray, fRhoValue, fhJetPt[fTriggerType][0], fhJetArea[fTriggerType][0], fhJetQA[fTriggerType][0]);
 671       if(fIsEmbedding)
 672         {
 673           if(fRunDLHJet) RunJetQA(fDLJetArray, 0, fhJetPt[fTriggerType][1], fhJetArea[fTriggerType][1], fhJetQA[fTriggerType][1]);
 674           if(fRunPLHJet)  RunJetQA(fPLJetArray, 0, fhJetPt[fTriggerType][2], fhJetArea[fTriggerType][2], fhJetQA[fTriggerType][2]);
 675         }
 676     }
 677
 678   if(!fIsEmbedding)
 679     {
 680       if(fRunTrkQA) RunTrackQA();
 681       if(fRunLeadTrkQA) RunLeadTrkQA();
 682     }
 683
 684   if(fStudyKtEffects) StudyKtEffects();
 685
 686   PostData(1, fOutputList);
 687   return;
 688 }
 689
 690 //________________________________________________________________________
 691 Int_t AliAnalysisTaskHJetDphi::FindSingleIncTrigger(const TClonesArray *trackArray, Double_t &trigPt, Double_t &trigPhi, Double_t &trigEta, const Int_t arrayType)
 692 {
 693   Int_t trigIndex = -1;
 694   trigPt = -1;
 695   trigPhi = -999;
 696   trigEta = -999;
 697
 698   // arrayType: 0 -> data, 1 -> embedding, 2 -> MC
 699   if(!trackArray) return 0;
 700
 701   Int_t nTT = 0, counter = 0;
 702   Int_t ntracks = trackArray->GetEntries();
 703   TArrayI arr;
 704   arr.Set(ntracks);
 705   for(Int_t it=0; it<ntracks; it++)
 706     {
 707       AliVParticle *t = dynamic_cast<AliVParticle*>(trackArray->At(it));
 708       if(!t || t->Charge()==0 || !AcceptTrack(t)) continue;
 709
 710       if(fAnaType==0 && arrayType==1 && t->GetLabel()==0) continue;
 711
 712       if(fAnaType==1 && arrayType<2 && !IsGoodAODtrack(t) ) continue;
 713
 714       Double_t pt = t->Pt();
 715       if(fTTtype==0)
 716         {
 717           if (pt<fTTMaxPt && pt>=fTTMinPt)
 718             {
 719               nTT++;
 720               arr.AddAt(it,counter);
 721               counter++;
 722             }
 723         }
 724     }
 725   arr.Set(counter);
 726
 727   if(fTTtype==0)
 728     {
 729       if(counter==0) trigIndex = -1;
 730       else if(counter==1) trigIndex = arr.At(0);
 731       else
 732         {
 733           fRandom->SetSeed(arr.At(0)); //make this random selection reproducible
 734           Double_t pro = fRandom->Uniform() * counter;
 735           trigIndex = arr.At(TMath::FloorNint(pro));
 736         }
 737       arr.Reset();
 738     }
 739
 740   if(trigIndex>-1)
 741     {
 742       AliVParticle *tt = (AliVParticle*) fTrackArray->At(trigIndex);
 743       trigPt = tt->Pt();
 744       trigPhi = tt->Phi();
 745       trigEta = tt->Eta();
 746
 747       if(fSwitchOnAvoidTpcHole)
 748         {
 749           if(fAvoidTpcHole==1 && !(trigPhi>3.89 && trigPhi<5.53)) trigIndex = -1;
 750           if(fAvoidTpcHole==2 && !(trigPhi>2.45 && trigPhi<3.44)) trigIndex = -1;
 751         }
 752
 753       if(fCutTPCBoundary)
 754         {
 755           Double_t phiDist = trigPhi - TMath::FloorNint(trigPhi/kSector)*kSector;
 756           if(phiDist<fDistToTPCBoundary || phiDist>kSector-fDistToTPCBoundary)
 757             {
 758               trigIndex = -1;
 759             }
 760         }
 761     }
 762
 763   if(trigIndex==-1) { trigPt = -1; trigPhi = -999; trigEta = -999;}
 764   if(arrayType<2) fTriggerTrkIndex = trigIndex;
 765   return nTT;
 766 }
 767
 768 //________________________________________________________________________
 769 void AliAnalysisTaskHJetDphi::RunSingleInclHJetCorr(Double_t trigPt, Double_t trigPhi, Double_t trigEta, const TClonesArray *jetArray, Double_t rho, THnSparse *hTT, THnSparse *hn)
 770 {
 771   if(trigPt<0 || !fJetArray) return;
 772
 773   if(hTT)
 774     {
 775       Double_t fillTT[] = {trigPt, fCentrality, (Double_t)fPtHardBin};
 776       hTT->Fill(fillTT);
 777     }
 778
 779   Int_t nJets = jetArray->GetEntries();
 780   Double_t jetPt = 0, jetEta = 0, jetPhi = 0, jetArea = 0;
 781   for(Int_t ij=0; ij<nJets; ij++)
 782     {
 783       if(fAnaType==0)
 784         {
 785           AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(jetArray->At(ij));
 786           jetPt   = jet->Pt();
 787           jetEta  = jet->Eta();
 788           jetPhi  = jet->Phi();
 789           jetArea = jet->Area();
 790         }
 791       else if(fAnaType==1)
 792         {
 793           AliAODJet* jet = dynamic_cast<AliAODJet*>(jetArray->At(ij));
 794           jetPt   = jet->Pt();
 795           jetEta  = jet->Eta();
 796           jetPhi  = jet->Phi();
 797           jetArea = jet->EffectiveAreaCharged();
 798         }
 799       else
 800         return;
 801       if(!IsGoodJet(jetEta)) continue; // eta cut
 802       Double_t dPhi = CalculateDPhi(trigPhi,jetPhi);
 803
 804       Double_t jetPtCorr = jetPt-rho*jetArea;
 805       if(jetPtCorr>fJetPtMin)
 806         {
 807           Double_t fill[] = {trigPt,jetPtCorr,dPhi,jetArea,fCentrality,trigEta-jetEta, (Double_t)fPtHardBin,static_cast<Double_t>(Entry()%10)};
 808           hn->Fill(fill);
 809         }
 810     }
 811 }
 812
 813
 814 //________________________________________________________________________
 815 void AliAnalysisTaskHJetDphi::RunTrackQA()
 816 {
 817   if(fIsEmbedding) return;
 818   if(!fTrackArray) return;
 819
 820   Int_t ntracks = fTrackArray->GetEntries();
 821   for(Int_t it=0; it<ntracks; it++)
 822     {
 823       AliVParticle *t = dynamic_cast<AliVParticle*>(fTrackArray->At(it));
 824       if(!t || t->Charge()==0 || !AcceptTrack(t)) continue;
 825       if(fAnaType==1 && !IsGoodAODtrack(t)) continue;
 826       fhTrkPt[fTriggerType]->Fill(t->Pt(),fCentrality);
 827       Double_t phi = t->Phi();
 828       if(phi<0) phi += 2*pi;
 829       Double_t fill[] = {t->Pt(), phi*TMath::RadToDeg(), t->Eta(), fCentrality};
 830       fhTrkQA[fTriggerType]->Fill(fill);
 831     }
 832
 833   if(fESD)
 834     {
 835       Int_t ntrack = fESD->GetNumberOfTracks();
 836       for(Int_t itr=0; itr<ntrack; itr++)
 837         {
 838           AliESDtrack *esdtrack = fESD->GetTrack(itr);
 839           if(!esdtrack ||  TMath::Abs(esdtrack->Eta())>0.9 )continue;
 840           Int_t type = -1;
 841           if(fEsdTrkCut->AcceptTrack(esdtrack))
 842             type = 0;
 843           else if(fEsdHybCut->AcceptTrack(esdtrack))
 844             type = 1;
 845           else
 846             continue;
 847
 848           Double_t resolution = -1;
 849           Int_t label = esdtrack->GetLabel();
 850           if(label>0 && fMC)
 851             {
 852               AliMCParticle *part = (AliMCParticle*)fMC->GetTrack(label);
 853               TParticle *tpart = (TParticle*)part->Particle();
 854               if(TMath::Abs(tpart->GetPdgCode())==211)
 855                 {
 856                   Double_t fillPhiRes[] = {esdtrack->Pt(),part->Phi()-esdtrack->Phi(),(Double_t)type,fCentrality};
 857                   fhTrkPhiRes[fTriggerType]->Fill(fillPhiRes);
 858                 }
 859               resolution = (part->Pt()-esdtrack->Pt())/part->Pt();
 860             }
 861           Double_t fillPtRes[] = {esdtrack->Pt(),esdtrack->Pt()*TMath::Sqrt(esdtrack->GetSigma1Pt2()),resolution,(Double_t)type,fCentrality};
 862           fhTrkPtRes[fTriggerType]->Fill(fillPtRes);
 863         }
 864     }
 865   else if(fAODIn)
 866     {
 867       ntracks = fAODIn->GetNumberOfTracks();
 868       Int_t type = -1;
 869       for(Int_t itrack=0; itrack<ntracks; itrack++)
 870         {
 871           AliAODTrack *aodtrack = (AliAODTrack*)fAODIn->GetTrack(itrack);
 872           if(!aodtrack || !AcceptTrack(dynamic_cast<AliVParticle*>(aodtrack)) ) continue;
 873           if (fAODfilterBits[0] < 0)
 874             {
 875               if (aodtrack->IsHybridGlobalConstrainedGlobal())
 876                 type = 3;
 877               else
 878                 continue;
 879             }
 880           else
 881             {
 882               if (aodtrack->TestFilterBit(fAODfilterBits[0]))  type = 0;
 883               else if (aodtrack->TestFilterBit(fAODfilterBits[1]) && (aodtrack->GetStatus()&AliESDtrack::kITSrefit)!=0) type = 1;
 884               else continue;
 885             }
 886           Double_t sigma1Pt2 = GetAODTrackPtRes(aodtrack);
 887           Double_t resolution = -1;
 888           Int_t label = aodtrack->GetLabel();
 889           if(label>0 && fMC)
 890             {
 891               AliMCParticle *part = (AliMCParticle*)fMC->GetTrack(label);
 892               resolution = (part->Pt()-aodtrack->Pt())/part->Pt();
 893               Double_t fillPhiRes[] = {aodtrack->Pt(),part->Phi()-aodtrack->Phi(),(Double_t)type,fCentrality};
 894               fhTrkPhiRes[fTriggerType]->Fill(fillPhiRes);
 895             }
 896           if(sigma1Pt2>0)
 897             {
 898               Double_t fillPtRes[5] = {aodtrack->Pt(),aodtrack->Pt()*TMath::Sqrt(sigma1Pt2),resolution,(Double_t)type,fCentrality};
 899               fhTrkPtRes[fTriggerType]->Fill(fillPtRes);
 900             }
 901         }
 902     }
 903 }
 904
 905 //________________________________________________________________________
 906 void AliAnalysisTaskHJetDphi::RunJetQA(const TClonesArray *jetArray, const Double_t rho, THnSparse *hJetPt, THnSparse *hJetArea, THnSparse *hJetQA)
 907 {
 908   Int_t nJets = jetArray->GetEntries();
 909   Double_t jetPt = 0, jetEta = 0, jetPhi = 0, jetArea = 0, jetPtCorr = 0;
 910   Int_t nCons = 0;
 911   for(Int_t ij=0; ij<nJets; ij++)
 912     {
 913       if(fAnaType==0)
 914         {
 915           AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(jetArray->At(ij));
 916           jetPt   = jet->Pt();
 917           jetEta  = jet->Eta();
 918           jetPhi  = jet->Phi();
 919           jetArea = jet->Area();
 920           nCons   = jet->GetNumberOfConstituents();
 921         }
 922       else if(fAnaType==1)
 923         {
 924           AliAODJet* jet = dynamic_cast<AliAODJet*>(jetArray->At(ij));
 925           jetPt   = jet->Pt();
 926           jetEta  = jet->Eta();
 927           jetPhi  = jet->Phi();
 928           jetArea = jet->EffectiveAreaCharged();
 929           nCons   = jet->GetRefTracks()->GetEntriesFast();
 930         }
 931       else
 932         return;
 933       if(!IsGoodJet(jetEta)) continue; // eta cut
 934
 935       jetPtCorr = jetPt-rho*jetArea;
 936       Double_t fillPt[] = {jetPtCorr, jetPt, fCentrality, (Double_t)fPtHardBin};
 937       hJetPt->Fill(fillPt);
 938
 939       Double_t fillA[] = {jetPtCorr, jetArea, fCentrality, (Double_t)fPtHardBin};
 940       hJetArea->Fill(fillA);
 941
 942       Double_t fillQA[] = {jetPtCorr, jetPhi*TMath::RadToDeg(), jetEta, jetArea, (Double_t)nCons, fCentrality, (Double_t)fPtHardBin};
 943       hJetQA->Fill(fillQA);
 944     }
 945 }
 946
 947
 948 //________________________________________________________________________
 949 void AliAnalysisTaskHJetDphi::RunLeadTrkQA()
 950 {
 951   if(fIsEmbedding || fTriggerTrkIndex<0) return;
 952   Double_t jetPt = -1, jetPhi = -999;
 953
 954   if(fAnaType==0)
 955     {
 956       Int_t nJets = fJetArray->GetEntries();
 957       for(Int_t ij=0; ij<nJets; ij++)
 958         {
 959           AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(fJetArray->At(ij));
 960           if(!IsGoodJet(jet->Eta())) continue; // eta cut
 961           for (Int_t i = 0; i < jet->GetNumberOfTracks(); i++)
 962             {
 963               if(jet->TrackAt(i)==fTriggerTrkIndex)
 964                 {
 965                   jetPt = jet->Pt();
 966                   jetPhi = jet->Phi();
 967                   break;
 968                 }
 969             }
 970         }
 971     }
 972   if(fAnaType==1)
 973     {
 974       jetPt = -1;
 975     }
 976
 977
 978   if(jetPt<=0) return;
 979   AliVParticle *tt = (AliVParticle*) fTrackArray->At(fTriggerTrkIndex);
 980   Double_t fill[] = {tt->Pt(), jetPt, tt->Phi()-jetPhi, tt->Pt()/jetPt, fCentrality};
 981   fhLeadTrkQA[fTriggerType]->Fill(fill);
 982 }
 983
 984
 985 //________________________________________________________________________
 986 void AliAnalysisTaskHJetDphi::StudyKtEffects()
 987 {
 988   if(fAnaType==1) return;
 989   if(fTriggerTrkIndex<0) return;
 990
 991   fKtValue = 0; // dummy
 992
 993
 994   AliPicoTrack *tt = (AliPicoTrack*) fTrackArray->At(fTriggerTrkIndex);
 995   Double_t triggerPhi = tt->GetTrackPhiOnEMCal();
 996
 997   Int_t nJets = fDLJetArray->GetEntries();
 998   for(Int_t ij=0; ij<nJets; ij++)
 999     {
1000       AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(fDLJetArray->At(ij));
1001       if(!IsGoodJet(jet->Eta())) continue; // eta cut
1002       Double_t jetPhi = jet->Phi();
1003       Double_t jetPt  = jet->Pt();
1004       Double_t dPhi = CalculateDPhi(triggerPhi,jetPhi);
1005       if(dPhi<pi+0.6 && dPhi>pi-0.6)
1006         {
1007           Double_t fill[] = {tt->Pt(), jetPt, jetPt*TMath::Tan(tt->GetTrackPtOnEMCal()), fCentrality, (Double_t)fPtHardBin};
1008           fhKtEffects[fTriggerType]->Fill(fill);
1009         }
1010     }
1011 }
1012
1013
1014 //________________________________________________________________________
1015 Bool_t AliAnalysisTaskHJetDphi::RetrieveArraies()
1016 {
1017   if(fAnaType==0)
1018     {
1019       // Get mc particles
1020       if (!fMcParticleArrName.IsNull())
1021         {
1022           AliInfo(Form("Retrieve mc particles %s!", fMcParticleArrName.Data()));
1023           fMcParticleArray = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fMcParticleArrName));
1024           if (!fMcParticleArray)
1025             {
1026               AliError(Form("Could not retrieve mc particles %s!", fMcParticleArrName.Data()));
1027               return kFALSE;
1028             }
1029
1030           TString fMcParticleMapArrName = fMcParticleArrName + "_Map";
1031           fMcParticlelMap = dynamic_cast<AliNamedArrayI*>(fEvent->FindListObject(fMcParticleMapArrName));
1032           if (!fMcParticlelMap)
1033             {
1034               AliWarning(Form("%s: Could not retrieve map for MC particles %s! Will assume MC labels consistent with indexes...", GetName(), fMcParticleArrName.Data()));
1035               return kFALSE;
1036             }
1037         }
1038
1039       // Get track collection
1040       if (!fTrackArrName.IsNull())
1041         {
1042           AliInfo(Form("Retrieve tracks %s!", fTrackArrName.Data()));
1043           fTrackArray = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fTrackArrName));
1044           if (!fTrackArray)
1045             {
1046               AliError(Form("Could not retrieve tracks %s!", fTrackArrName.Data()));
1047               return kFALSE;
1048             }
1049         }
1050
1051       // Get mixed track collection
1052       if (!fEmbTrkArrName.IsNull())
1053         {
1054           AliInfo(Form("Retrieve PYTHIA+PbPb tracks %s!", fEmbTrkArrName.Data()));
1055           fEmbTrkArray = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fEmbTrkArrName));
1056           if (!fEmbTrkArray)
1057             {
1058               AliError(Form("Could not retrieve PYTHIA+PbPb tracks %s!", fEmbTrkArrName.Data()));
1059               return kFALSE;
1060             }
1061         }
1062
1063       // Get Rho value
1064       if(fCollisionSystem=="PbPb")
1065         {
1066           if(!fRhoName.IsNull())
1067             {
1068               AliInfo(Form("Retrieve rho %s!", fRhoName.Data()));
1069               fRho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRhoName));
1070               if(!fRho)
1071                 {
1072                   AliError(Form("Could not retrieve rho %s!",fRhoName.Data()));
1073                   return kFALSE;
1074                 }
1075             }
1076         }
1077
1078       // Get jet collection
1079       if (!fJetArrName.IsNull())
1080         {
1081           AliInfo(Form("Retrieve jets %s!", fJetArrName.Data()));
1082           fJetArray = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fJetArrName));
1083           if (!fJetArray)
1084             {
1085               AliError(Form("%s: Could not retrieve jets %s!", GetName(), fJetArrName.Data()));
1086               return kFALSE;
1087             }
1088         }
1089
1090       // Get DL jet collection
1091       if (!fDLJetArrName.IsNull())
1092         {
1093           AliInfo(Form("Retrieve DL jets %s!", fDLJetArrName.Data()));
1094           fDLJetArray = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fDLJetArrName));
1095           if (!fDLJetArray)
1096             {
1097               AliError(Form("%s: Could not retrieve DL jets %s!", GetName(), fDLJetArrName.Data()));
1098               return kFALSE;
1099             }
1100         }
1101
1102       // Get PL jet collection
1103       if (!fPLJetArrName.IsNull())
1104         {
1105           AliInfo(Form("Retrieve PL jets %s!", fPLJetArrName.Data()));
1106           fPLJetArray = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fPLJetArrName));
1107           if (!fPLJetArray)
1108             {
1109               AliError(Form("%s: Could not retrieve PL jets %s!", GetName(), fPLJetArrName.Data()));
1110               return kFALSE;
1111             }
1112         }
1113
1114       if(fIsEmbedding)
1115         {
1116           if(fAnaType==0 && !fPtHardBinName)
1117             {
1118               // Get embedded pt hard bin number
1119               fPtHardBinName = static_cast<AliNamedString*>(fEvent->FindListObject("AODEmbeddingFile"));
1120               if(!fPtHardBinName)
1121                 {
1122                   AliError("The object for pt hard bin information is not available!");
1123                   return kFALSE;
1124                 }
1125             }
1126         }
1127     }
1128
1129   if(fAnaType==1)
1130     {
1131       // Get mc particles
1132       if (!fMcParticleArrName.IsNull())
1133         {
1134           AliInfo(Form("Retrieve mc particles %s!", fMcParticleArrName.Data()));
1135           if(fAODOut        && !fMcParticleArray) fMcParticleArray = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fMcParticleArrName));
1136           if(fAODIn         && !fMcParticleArray) fMcParticleArray = dynamic_cast<TClonesArray*>(fAODIn ->FindListObject(fMcParticleArrName));
1137           if(fAODExtension  && !fMcParticleArray) fMcParticleArray = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fMcParticleArrName));
1138           if (!fMcParticleArray)
1139             {
1140               AliError(Form("Could not retrieve mc particles %s!", fMcParticleArrName.Data()));
1141               return kFALSE;
1142             }
1143         }
1144
1145       // Get tracks
1146       if (!fTrackArrName.IsNull())
1147         {
1148           AliInfo(Form("Retrieve tracks %s!", fTrackArrName.Data()));
1149           if(fAODIn         && !fTrackArray) fTrackArray = dynamic_cast<TClonesArray*>(fAODIn ->FindListObject(fTrackArrName));
1150           if(fAODOut        && !fTrackArray) fTrackArray = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fTrackArrName));
1151           if (!fTrackArray)
1152             {
1153               AliError(Form("Could not retrieve tracks %s!", fTrackArrName.Data()));
1154               return kFALSE;
1155             }
1156         }
1157
1158       // Get PYTHIA+PbPb tracks
1159       if (!fEmbTrkArrName.IsNull())
1160         {
1161           AliInfo(Form("Retrieve PYTHIA+PbPb tracks %s!", fEmbTrkArrName.Data()));
1162           if(fAODOut        && !fEmbTrkArray) fEmbTrkArray = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fEmbTrkArrName));
1163           if(fAODIn         && !fEmbTrkArray) fEmbTrkArray = dynamic_cast<TClonesArray*>(fAODIn ->FindListObject(fEmbTrkArrName));
1164           if(fAODExtension  && !fEmbTrkArray) fEmbTrkArray = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fEmbTrkArrName));
1165           if (!fTrackArray)
1166             {
1167               AliError(Form("Could not retrieve PYTHIA+PbPb tracks %s!", fTrackArrName.Data()));
1168               return kFALSE;
1169             }
1170         }
1171
1172       // Get jet collection
1173       if (!fJetArrName.IsNull())
1174         {
1175           AliInfo(Form("Retrieve jets %s!", fJetArrName.Data()));
1176           if(fAODOut        && !fJetArray) fJetArray = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetArrName));
1177           if(fAODExtension  && !fJetArray) fJetArray = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetArrName));
1178           if(fAODIn         && !fJetArray) fJetArray = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetArrName));
1179           if (!fJetArray)
1180             {
1181               AliError(Form("%s: Could not retrieve jets %s!", GetName(), fJetArrName.Data()));
1182               return kFALSE;
1183             }
1184         }
1185
1186       // Get DL jet collection
1187       if (!fDLJetArrName.IsNull())
1188         {
1189           AliInfo(Form("Retrieve DL jets %s!", fDLJetArrName.Data()));
1190           if(fAODOut        && !fDLJetArray) fDLJetArray = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fDLJetArrName));
1191           if(fAODExtension  && !fDLJetArray) fDLJetArray = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fDLJetArrName));
1192           if(fAODIn         && !fDLJetArray) fDLJetArray = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fDLJetArrName));
1193           if (!fDLJetArray)
1194             {
1195               AliError(Form("%s: Could not retrieve DL jets %s!", GetName(), fDLJetArrName.Data()));
1196               return kFALSE;
1197             }
1198         }
1199
1200       // Get PL jet collection
1201       if (!fPLJetArrName.IsNull())
1202         {
1203           AliInfo(Form("Retrieve PL jets %s!", fPLJetArrName.Data()));
1204           if(fAODOut        && !fPLJetArray) fPLJetArray = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fPLJetArrName));
1205           if(fAODExtension  && !fPLJetArray) fPLJetArray = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fPLJetArrName));
1206           if(fAODIn         && !fPLJetArray) fPLJetArray = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fPLJetArrName));
1207           if (!fPLJetArray)
1208             {
1209               AliError(Form("%s: Could not retrieve PL jets %s!", GetName(), fPLJetArrName.Data()));
1210               return kFALSE;
1211             }
1212         }
1213
1214       // Get Rho
1215       if(fCollisionSystem=="PbPb" && !fRhoName.IsNull())
1216         {
1217           AliInfo(Form("Retrieve rho %s!", fRhoName.Data()));
1218           if(fAODOut       && !fEvtBkg ) fEvtBkg = (AliAODJetEventBackground*)(fAODOut->FindListObject(fRhoName));
1219           if(fAODExtension && !fEvtBkg ) fEvtBkg = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fRhoName));
1220           if(fAODIn        && !fEvtBkg ) fEvtBkg = (AliAODJetEventBackground*)(fAODIn->FindListObject(fRhoName));
1221           if(!fEvtBkg)
1222             {
1223               AliError(Form("Could not retrieve rho %s!",fRhoName.Data()));
1224               return kFALSE;
1225             }
1226         }
1227     }
1228
1229   return kTRUE;
1230 }
1231
1232 //________________________________________________________________________
1233 Double_t AliAnalysisTaskHJetDphi::CalculatePhi(const Double_t py, const Double_t px)
1234 {
1235   Double_t phi = TMath::ATan2(py,px);
1236   if(phi<0) phi += 2*pi;
1237   return phi;
1238 }
1239
1240 //________________________________________________________________________
1241 Double_t AliAnalysisTaskHJetDphi::CalculateDPhi(const Double_t phi1, const Double_t phi2)
1242 {
1243   Double_t dPhi = phi1-phi2;
1244   if(dPhi>2*pi)  dPhi -= 2*pi;
1245   if(dPhi<-2*pi) dPhi += 2*pi;
1246   if(dPhi<-0.5*pi) dPhi += 2*pi;
1247   if(dPhi>1.5*pi)  dPhi -= 2*pi;
1248   return dPhi;
1249 }
1250
1251 //________________________________________________________________________
1252 Int_t AliAnalysisTaskHJetDphi::LocateToTPCHole(const Double_t phi)
1253 {
1254   if(phi<4-pi/2 && phi>5.5-1.5*pi) return 1; // away-side
1255   else return 0;   //near-side
1256 }
1257
1258
1259 //________________________________________________________________________
1260 Bool_t AliAnalysisTaskHJetDphi::AcceptTrack(AliVParticle *track)
1261 {
1262   if(track->Pt()<fMinTrkPt || track->Pt()>fMaxTrkPt) return kFALSE;
1263   if(track->Eta()<fMinTrkEta || track->Eta()>fMaxTrkEta) return kFALSE;
1264   if(track->Phi()<fMinTrkPhi || track->Phi()>fMaxTrkPhi) return kFALSE;
1265   return kTRUE;
1266 }
1267
1268 //________________________________________________________________________
1269 Bool_t AliAnalysisTaskHJetDphi::IsGoodAODtrack(AliVParticle *track)
1270 {
1271   AliAODTrack *aodtrack = static_cast<AliAODTrack*>(track);
1272   if( fFilterMask>0)
1273     {
1274       if(!aodtrack->TestFilterBit(fFilterMask) ) return kFALSE;
1275     }
1276   else
1277     {
1278       if(!aodtrack->IsHybridGlobalConstrainedGlobal()) return kFALSE;
1279     }
1280   if( fRequireITSRefit && (aodtrack->GetStatus()&AliESDtrack::kITSrefit)==0 ) return kFALSE;
1281   if (fRequireSharedClsCut)
1282     {
1283       Double_t frac = Double_t(aodtrack->GetTPCnclsS())/Double_t(aodtrack->GetTPCncls());
1284       if (frac > 0.4) return kFALSE;
1285     }
1286   return kTRUE;
1287 }
1288
1289 //________________________________________________________________________
1290 Bool_t AliAnalysisTaskHJetDphi::IsGoodJet(Double_t jetEta)
1291 {
1292   Double_t etaCut = (0.9-fRadius>0.5)?0.5:0.9-fRadius;
1293   if(TMath::Abs(jetEta)>etaCut) return kFALSE;
1294   return kTRUE;
1295 }
1296
1297 //________________________________________________________________________
1298 Int_t AliAnalysisTaskHJetDphi::GetPtHardBin(Double_t ptHard)
1299 {
1300   const Int_t nBins = 10;
1301   Double_t binLimits[nBins] = { 5., 11., 21., 36., 57., 84., 117., 156., 200., 249. }; // lower limits
1302   Int_t bin = -1;
1303   while(bin<nBins-1 && binLimits[bin+1]<ptHard)
1304     {
1305       bin++;
1306     }
1307   return bin;
1308 }
1309
1310 //________________________________________________________________________
1311 Int_t AliAnalysisTaskHJetDphi::GetParticleType(Int_t pdg_input)
1312 {
1313   Int_t type = -1;
1314   Int_t pdg = TMath::Abs(pdg_input);
1315   if(pdg==211) type = 0;
1316   else if(pdg==321) type = 1;
1317   else if(pdg==2212) type = 2;
1318   else if(pdg==11) type = 3;
1319   else if(pdg>=3122 && pdg<=3334) type = 4;
1320   else type = 9;
1321   return type;
1322 }
1323
1324 //________________________________________________________________________
1325 Double_t AliAnalysisTaskHJetDphi::GetAODTrackPtRes(AliAODTrack *aodtrack)
1326 {
1327   Double_t sigma1Pt2 = -1;
1328   Double_t cov[21] = {0,}, pxpypz[3] = {0,}, xyz[3] = {0,};
1329   AliExternalTrackParam *exParam = new  AliExternalTrackParam();
1330   aodtrack->GetCovMatrix(cov);
1331   aodtrack->PxPyPz(pxpypz);
1332   aodtrack->GetXYZ(xyz);
1333   exParam->Set(xyz,pxpypz,cov,aodtrack->Charge());
1334   sigma1Pt2 = exParam->GetSigma1Pt2();
1335   delete exParam;
1336   return sigma1Pt2;
1337 }
1338
1339 //
1340 //________________________________________________________________________
1341 //
1342 void AliAnalysisTaskHJetDphi::PrintConfig()
1343 {
1344   const char *decision[2] = {"no","yes"};
1345   printf("\n\n===== h-jet analysis configuration =====\n");
1346   printf("Input event type: %s - %s\n",fCollisionSystem.Data(),fPeriod.Data());
1347   printf("Is this MC data: %s\n",decision[fIsMC]);
1348   printf("Run on particle level: %s\n",decision[fAnalyzeMCTruth]);
1349   printf("Event type selection: %d\n",fOfflineTrgMask);
1350   printf("Is embedding? %s\n",decision[fIsEmbedding]);
1351   printf("Track filter mask: %d\n",fFilterMask);
1352   printf("Require track to have ITS refit? %s\n",decision[fRequireITSRefit]);
1353   printf("Require to cut on fraction of shared TPC clusters? %s\n",decision[fRequireSharedClsCut]);
1354   printf("Track pt range: %2.2f < pt < %2.2f\n",fMinTrkPt, fMaxTrkPt);
1355   printf("Track eta range: %2.1f < eta < %2.1f\n",fMinTrkEta, fMaxTrkEta);
1356   printf("Track phi range: %2.0f < phi < %2.0f\n",fMinTrkPhi*TMath::RadToDeg(),fMaxTrkPhi*TMath::RadToDeg());
1357   printf("Cut TT away from boundary: %s with distance = %2.2f\n",decision[fCutTPCBoundary],fDistToTPCBoundary);
1358   printf("Avoid TPC holes: %s\n", decision[fSwitchOnAvoidTpcHole]);
1359   printf("Jet cone size R = %1.1f, and jet pt > %1.0f GeV/c \n",fRadius,fJetPtMin);
1360   printf("Run track QA: %s\n",decision[fRunTrkQA]);
1361   printf("Run jet QA: %s\n",decision[fRunJetQA]);
1362   printf("Run single inclusive h+jet analysis: %s\n",decision[fRunSingleInclHJet]);
1363   printf("TT interval:    %2.0f < pt < %2.0f\n",fTTMinPt, fTTMaxPt);
1364   printf("Run leading track QA: %s\n",decision[fRunLeadTrkQA]);
1365   printf("Run kT effects study: %s\n",decision[fStudyKtEffects]);
1366   printf("Run background flow: %s\n",decision[fRunBkgFlow]);
1367   printf("=======================================\n\n");
1368 }
1369
1370 //________________________________________________________________________
1371 Double_t AliAnalysisTaskHJetDphi::GetZ(const Double_t trkPx, const Double_t trkPy, const Double_t trkPz, const Double_t jetPx, const Double_t jetPy, const Double_t jetPz)
1372 {
1373   return (trkPx*jetPx+trkPy*jetPy+trkPz*jetPz)/(jetPx*jetPx+jetPy*jetPy+jetPz*jetPz);
1374 }
1375
1376 //________________________________________________________________________
1377 void AliAnalysisTaskHJetDphi::Terminate(Option_t *)
1378 {
1379   // Called once at the end of the query
1380 }