PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetHadEPpid.cxx

   1 // Jet-Hadron Correlations PID
   2 // Event plane dependence task.
   3 //
   4 // Author: J Mazer
   5
   6 // task head include
   7 #include "AliAnalysisTaskEmcalJetHadEPpid.h"
   8
   9 // general ROOT includes
  10 #include <TCanvas.h>
  11 #include <TChain.h>
  12 #include <TClonesArray.h>
  13 #include <TH1F.h>
  14 #include <TH2F.h>
  15 #include <TH3F.h>
  16 #include <THnSparse.h>
  17 #include <TList.h>
  18 #include <TLorentzVector.h>
  19 #include <TParameter.h>
  20 #include <TParticle.h>
  21 #include <TTree.h>
  22 #include <TVector3.h>
  23 #include <TObjArray.h>
  24
  25 // AliROOT includes
  26 #include "AliAODEvent.h"
  27 #include "AliESDEvent.h"
  28 #include "AliAnalysisManager.h"
  29 #include "AliAnalysisTask.h"
  30 #include "AliCentrality.h"
  31 #include "AliEmcalJet.h"
  32 #include "AliAODJet.h"
  33 #include "AliVCluster.h"
  34 #include "AliVTrack.h"
  35 #include <AliVEvent.h>
  36 #include <AliVParticle.h>
  37 #include "AliRhoParameter.h"
  38 #include "AliEmcalParticle.h"
  39
  40 // Localized Rho includes
  41 #include "AliLocalRhoParameter.h"
  42 #include "AliAnalysisTaskLocalRho.h"
  43
  44 // event handler (and pico's) includes
  45 #include <AliInputEventHandler.h>
  46 #include <AliVEventHandler.h>
  47 #include "AliESDInputHandler.h"
  48 #include "AliPicoTrack.h"
  49 #include "AliEventPoolManager.h"
  50
  51 // PID includes
  52 #include "AliPIDResponse.h"
  53 #include "AliTPCPIDResponse.h"
  54 #include "AliESDpid.h"
  55
  56 // magnetic field includes
  57 #include "TGeoGlobalMagField.h"
  58 #include "AliMagF.h"
  59
  60 ClassImp(AliAnalysisTaskEmcalJetHadEPpid)
  61
  62 //________________________________________________________________________
  63 AliAnalysisTaskEmcalJetHadEPpid::AliAnalysisTaskEmcalJetHadEPpid() :
  64   AliAnalysisTaskEmcalJet("correlations",kFALSE),
  65   fPhimin(-10), fPhimax(10),
  66   fEtamin(-0.9), fEtamax(0.9),
  67   fAreacut(0.0), fTrkBias(5), fClusBias(5), fTrkEta(0.9), fJetPtcut(15.0),
  68   fDoEventMixing(0), fMixingTracks(50000),
  69   doPlotGlobalRho(kFALSE), doVariableBinning(kFALSE),
  70   dovarbinTHnSparse(kFALSE), useAOD(kTRUE), doPID(kTRUE),
  71   fLocalRhoVal(0),
  72   fTracksName(""), fJetsName(""),
  73   event(0),
  74   isPItpc(0), isKtpc(0), isPtpc(0), nPIDtpc(0) ,   // pid TPC
  75   isPIits(0), isKits(0), isPits(0), nPIDits(0) ,  // pid ITS
  76   isPItof(0), isKtof(0), isPtof(0), nPIDtof(0) ,  // pid TOF
  77   fPoolMgr(0x0),
  78   fPIDResponse(0x0), fTPCResponse(),
  79   fESD(0), fAOD(0),
  80   fHistTPCdEdX(0), fHistITSsignal(0), //fHistTOFsignal(0),
  81   fHistRhovsCent(0), fHistNjetvsCent(0), fHistCentrality(0),
  82   fHistZvtx(0), fHistMult(0),
  83   fHistJetPhi(0), fHistTrackPhi(0),
  84   fHistJetHaddPhiIN(0), fHistJetHaddPhiOUT(0), fHistJetHaddPhiMID(0),
  85   fHistJetHaddPhiBias(0), fHistJetHaddPhiINBias(0), fHistJetHaddPhiOUTBias(0), fHistJetHaddPhiMIDBias(0),
  86   fHistMEdPHI(0),
  87   fHistTrackPtallcent(0),
  88   fHistJetEtaPhi(0), fHistJetHEtaPhi(0),
  89   fHistSEphieta(0), fHistMEphieta(0),
  90   fHistJetHaddPHI(0),
  91   fHistPID(0),
  92   fhnPID(0x0), fhnMixedEvents(0x0), fhnJH(0x0)
  93 {
  94   // Default Constructor
  95   for(Int_t ilab=0; ilab<4; ilab++){
  96     for(Int_t ipta=0; ipta<7; ipta++){
  97       fHistTrackEtaPhi[ilab][ipta]=0;
  98     }  // end of pt-associated loop
  99   } // end of lab loop
 100
 101   for(Int_t itrackpt=0; itrackpt<9; itrackpt++){
 102     fHistJetHadbindPhi[itrackpt]=0;
 103     fHistJetHadbindPhiIN[itrackpt]=0;
 104     fHistJetHadbindPhiMID[itrackpt]=0;
 105     fHistJetHadbindPhiOUT[itrackpt]=0;
 106   } // end of trackpt bin loop
 107
 108   for(Int_t icent = 0; icent<2; ++icent){
 109     fHistJetPt[icent]=0;
 110     fHistJetPtBias[icent]=0;
 111     fHistJetPtTT[icent]=0;
 112     fHistAreavsRawPt[icent]=0;
 113
 114     for(Int_t iptjet = 0; iptjet<5; ++iptjet){
 115       for(Int_t ieta = 0; ieta<3; ++ieta){
 116         fHistJetH[icent][iptjet][ieta]=0;
 117         fHistJetHBias[icent][iptjet][ieta]=0;
 118         fHistJetHTT[icent][iptjet][ieta]=0;
 119       } // end of eta loop
 120     } // end of pt-jet loop
 121   } // end of centrality loop
 122
 123   // centrality dependent histo's
 124   for (Int_t i = 0;i<6;++i){
 125     fHistJetPtvsTrackPt[i]      = 0;
 126     fHistRawJetPtvsTrackPt[i]   = 0;
 127     fHistTrackPt[i]             = 0;
 128     fHistEP0[i]                 = 0;
 129     fHistEP0A[i]                = 0;
 130     fHistEP0C[i]                = 0;
 131     fHistEPAvsC[i]              = 0;
 132     fHistJetPtvsdEP[i]          = 0;
 133     fHistJetPtvsdEPBias[i]      = 0;
 134     fHistRhovsdEP[i]            = 0;
 135     fHistJetEtaPhiPt[i]         = 0;
 136     fHistJetEtaPhiPtBias[i]     = 0;
 137     fHistJetPtArea[i]           = 0;
 138     fHistJetPtAreaBias[i]       = 0;
 139     fHistJetPtNcon[i]           = 0;
 140     fHistJetPtNconBias[i]       = 0;
 141     fHistJetPtNconCh[i]         = 0;
 142     fHistJetPtNconBiasCh[i]     = 0;
 143     fHistJetPtNconEm[i]         = 0;
 144     fHistJetPtNconBiasEm[i]     = 0;
 145     fHistJetHaddPhiINcent[i]    = 0;
 146     fHistJetHaddPhiOUTcent[i]   = 0;
 147     fHistJetHaddPhiMIDcent[i]   = 0;
 148   }
 149
 150   SetMakeGeneralHistograms(kTRUE);
 151
 152   // define input and output slots here
 153   DefineInput(0, TChain::Class());
 154   DefineOutput(1, TList::Class());
 155 }
 156
 157 //________________________________________________________________________
 158 AliAnalysisTaskEmcalJetHadEPpid::AliAnalysisTaskEmcalJetHadEPpid(const char *name) :
 159   AliAnalysisTaskEmcalJet(name,kTRUE),
 160   fPhimin(-10), fPhimax(10),
 161   fEtamin(-0.9), fEtamax(0.9),
 162   fAreacut(0.0), fTrkBias(5), fClusBias(5), fTrkEta(0.9), fJetPtcut(15.0),
 163   fDoEventMixing(0), fMixingTracks(50000),
 164   doPlotGlobalRho(kFALSE), doVariableBinning(kFALSE),
 165   dovarbinTHnSparse(kFALSE), useAOD(kTRUE), doPID(kTRUE),
 166   fLocalRhoVal(0),
 167   fTracksName(""), fJetsName(""),
 168   event(0),
 169   isPItpc(0), isKtpc(0), isPtpc(0), nPIDtpc(0),    // pid TPC
 170   isPIits(0), isKits(0), isPits(0), nPIDits(0),   // pid ITS
 171   isPItof(0), isKtof(0), isPtof(0), nPIDtof(0),   // pid TOF
 172   fPoolMgr(0x0),
 173   fPIDResponse(0x0), fTPCResponse(),
 174   fESD(0), fAOD(0),
 175   fHistTPCdEdX(0), fHistITSsignal(0), //fHistTOFsignal(0),
 176   fHistRhovsCent(0), fHistNjetvsCent(0), fHistCentrality(0),
 177   fHistZvtx(0), fHistMult(0),
 178   fHistJetPhi(0), fHistTrackPhi(0),
 179   fHistJetHaddPhiIN(0), fHistJetHaddPhiOUT(0), fHistJetHaddPhiMID(0),
 180   fHistJetHaddPhiBias(0), fHistJetHaddPhiINBias(0), fHistJetHaddPhiOUTBias(0), fHistJetHaddPhiMIDBias(0),
 181   fHistMEdPHI(0),
 182   fHistTrackPtallcent(0),
 183   fHistJetEtaPhi(0), fHistJetHEtaPhi(0),
 184   fHistSEphieta(0), fHistMEphieta(0),
 185   fHistJetHaddPHI(0),
 186   fHistPID(0),
 187   fhnPID(0x0), fhnMixedEvents(0x0), fhnJH(0x0)
 188 {
 189   // Default Constructor
 190   for(Int_t ilab=0; ilab<4; ilab++){
 191     for(Int_t ipta=0; ipta<7; ipta++){
 192       fHistTrackEtaPhi[ilab][ipta]=0;
 193     }  // end of pt-associated loop
 194   } // end of lab loop
 195
 196   for(Int_t itrackpt=0; itrackpt<9; itrackpt++){
 197     fHistJetHadbindPhi[itrackpt]=0;
 198     fHistJetHadbindPhiIN[itrackpt]=0;
 199     fHistJetHadbindPhiMID[itrackpt]=0;
 200     fHistJetHadbindPhiOUT[itrackpt]=0;
 201   } // end of trackpt bin loop
 202
 203   for(Int_t icent = 0; icent<2; ++icent){
 204     fHistJetPt[icent]=0;
 205     fHistJetPtBias[icent]=0;
 206     fHistJetPtTT[icent]=0;
 207     fHistAreavsRawPt[icent]=0;
 208
 209     for(Int_t iptjet = 0; iptjet<5; ++iptjet){
 210       for(Int_t ieta = 0; ieta<3; ++ieta){
 211         fHistJetH[icent][iptjet][ieta]=0;
 212         fHistJetHBias[icent][iptjet][ieta]=0;
 213         fHistJetHTT[icent][iptjet][ieta]=0;
 214       } // end of eta loop
 215     } // end of pt-jet loop
 216   } // end of centrality loop
 217
 218   // centrality dependent histo's
 219   for (Int_t i = 0;i<6;++i){
 220     fHistJetPtvsTrackPt[i]      = 0;
 221     fHistRawJetPtvsTrackPt[i]   = 0;
 222     fHistTrackPt[i]             = 0;
 223     fHistEP0[i]                 = 0;
 224     fHistEP0A[i]                = 0;
 225     fHistEP0C[i]                = 0;
 226     fHistEPAvsC[i]              = 0;
 227     fHistJetPtvsdEP[i]          = 0;
 228     fHistJetPtvsdEPBias[i]      = 0;
 229     fHistRhovsdEP[i]            = 0;
 230     fHistJetEtaPhiPt[i]         = 0;
 231     fHistJetEtaPhiPtBias[i]     = 0;
 232     fHistJetPtArea[i]           = 0;
 233     fHistJetPtAreaBias[i]       = 0;
 234     fHistJetPtNcon[i]           = 0;
 235     fHistJetPtNconBias[i]       = 0;
 236     fHistJetPtNconCh[i]         = 0;
 237     fHistJetPtNconBiasCh[i]     = 0;
 238     fHistJetPtNconEm[i]         = 0;
 239     fHistJetPtNconBiasEm[i]     = 0;
 240     fHistJetHaddPhiINcent[i]    = 0;
 241     fHistJetHaddPhiOUTcent[i]   = 0;
 242     fHistJetHaddPhiMIDcent[i]   = 0;
 243   }
 244
 245   SetMakeGeneralHistograms(kTRUE);
 246
 247   // define input and output slots here
 248   DefineInput(0, TChain::Class());
 249   DefineOutput(1, TList::Class());
 250 }
 251
 252 //_______________________________________________________________________
 253 AliAnalysisTaskEmcalJetHadEPpid::~AliAnalysisTaskEmcalJetHadEPpid()
 254 {
 255   // destructor
 256   if (fOutput) {
 257     delete fOutput;
 258     fOutput = 0;
 259   }
 260 }
 261
 262 //________________________________________________________________________
 263 void AliAnalysisTaskEmcalJetHadEPpid::UserCreateOutputObjects()
 264 { // This is called ONCE!
 265   if (!fCreateHisto) return;
 266   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 267   OpenFile(1); // do I need the 1?
 268
 269   // create histograms that arn't array
 270   fHistTPCdEdX = new TH2F("TPCdEdX", "TPCdEdX", 2000, 0.0, 100.0, 500, 0, 500);
 271   fHistITSsignal = new TH2F("ITSsignal", "ITSsignal", 2000, 0.0, 100.0, 500, 0, 500);
 272 //  fHistTOFsignal = new TH2F("TOFsignal", "TOFsignal", 2000, 0.0, 100.0, 500, 0, 500);
 273   fHistRhovsCent             = new TH2F("RhovsCent",              "RhovsCent",             100, 0.0, 100.0, 500, 0, 500);
 274   fHistNjetvsCent            = new TH2F("NjetvsCent",             "NjetvsCent",            100, 0.0, 100.0, 100, 0, 100);
 275   fHistTrackPtallcent = new TH1F("fHistTrackPtallcent", "p_{T} distribution", 1000, 0.0, 100.0);
 276   fHistCentrality = new TH1F("fHistCentrality","centrality",100,0,100);
 277   fHistZvtx = new TH1F("fHistZvertex","z vertex",60,-30,30);
 278   fHistJetEtaPhi = new TH2F("fHistJetEtaPhi","Jet #eta-#phi",900,-1.8,1.8,640,-3.2,3.2);
 279   fHistJetHEtaPhi = new TH2F("fHistJetHEtaPhi","Jet-Hadron #Delta#eta-#Delta#phi",900,-1.8,1.8,640,-1.6,4.8);
 280   fHistSEphieta = new TH2F("fHistSEphieta", "Single Event #phi-#eta distribution", 64,-0.5*TMath::Pi(), 1.5*TMath::Pi(), 32,-1.8,1.8);  // was 64 bins
 281   fHistMEphieta = new TH2F("fHistMEphieta", "Mixed Event #phi-#eta distribution", 64, -0.5*TMath::Pi(), 1.5*TMath::Pi(), 32,-1.8,1.8);  // was 64 bins
 282   fHistJetHaddPHI = new TH1F("fHistJetHaddPHI", "Jet-Hadron #Delta#varphi", 128,-0.5*TMath::Pi(),1.5*TMath::Pi());
 283   fHistJetPhi = new TH1F("fHistJetPhi", "Jet #phi Distribution", 128, -2.0*TMath::Pi(), 2.0*TMath::Pi());
 284   fHistTrackPhi = new TH1F("fHistTrackPhi", "Track #phi Distribution", 128, -2.0*TMath::Pi(), 2.0*TMath::Pi());
 285   fHistJetHaddPhiIN = new TH1F("fHistJetHaddPhiIN","Jet-Hadron #Delta#varphi IN PLANE", 128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
 286   fHistJetHaddPhiOUT = new TH1F("fHistJetHaddPhiOUT","Jet-Hadron #Delta#varphi OUT PLANE",128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
 287   fHistJetHaddPhiMID = new TH1F("fHistJetHaddPhiMID","Jet-Hadron #Delta#varphi MIDDLE of PLANE",128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
 288   fHistJetHaddPhiBias = new TH1F("fHistJetHaddPhiBias","Jet-Hadron #Delta#varphi with bias",128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
 289   fHistJetHaddPhiINBias = new TH1F("fHistJetHaddPhiINBias","Jet-Hadron #Delta#varphi IN PLANE with bias", 128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
 290   fHistJetHaddPhiOUTBias = new TH1F("fHistJetHaddPhiOUTBias","Jet-Hadron #Delta#varphi OUT PLANE with bias",128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
 291   fHistJetHaddPhiMIDBias = new TH1F("fHistJetHaddPhiMIDBias","Jet-Hadron #Delta#varphi MIDDLE of PLANE with bias",128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
 292
 293   fHistPID = new TH1F("fHistPID","PID Counter",12, -0.5, 11.5);
 294   SetfHistPIDcounterLabels(fHistPID);
 295   fOutput->Add(fHistPID);
 296
 297   char name[200];
 298
 299   for(Int_t ipta=0; ipta<7; ++ipta){
 300     for(Int_t ilab=0; ilab<4; ++ilab){
 301       snprintf(name,200, "fHistTrackEtaPhi_%i_%i", ilab,ipta);
 302       fHistTrackEtaPhi[ilab][ipta] = new TH2F(name,name,400,-1,1,640,0.0,2.*TMath::Pi());
 303 ///      fOutput->Add(fHistTrackEtaPhi[ilab][ipta]);
 304     } // end of lab loop
 305   } // end of pt-associated loop
 306
 307
 308   // jet hadron (binned) correlations in dPHI
 309   for(Int_t itrackpt=0; itrackpt<9; itrackpt++){
 310     snprintf(name,200,"fHistJetHadbindPhi_%i",itrackpt);
 311     fHistJetHadbindPhi[itrackpt] = new TH1F(name,name,128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
 312     fOutput->Add(fHistJetHadbindPhi[itrackpt]);
 313
 314     snprintf(name,200,"fHistJetHadbindPhiIN_%i",itrackpt);
 315     fHistJetHadbindPhiIN[itrackpt] = new TH1F(name,name,128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
 316     fOutput->Add(fHistJetHadbindPhiIN[itrackpt]);
 317
 318     snprintf(name,200,"fHistJetHadbindPhiMID_%i",itrackpt);
 319     fHistJetHadbindPhiMID[itrackpt] = new TH1F(name,name,128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
 320     fOutput->Add(fHistJetHadbindPhiMID[itrackpt]);
 321
 322     snprintf(name,200,"fHistJetHadbindPhiOUT_%i",itrackpt);
 323     fHistJetHadbindPhiOUT[itrackpt] = new TH1F(name,name,128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
 324     fOutput->Add(fHistJetHadbindPhiOUT[itrackpt]);
 325   } // end of trackpt bin loop
 326
 327   for(Int_t icent = 0; icent<2; ++icent){
 328     snprintf(name,200,"fHistJetPt_%i",icent);
 329     fHistJetPt[icent] = new TH1F(name,name,200,0,200);
 330     fOutput->Add(fHistJetPt[icent]);
 331
 332     snprintf(name,200,"fHistJetPtBias_%i",icent);
 333     fHistJetPtBias[icent] = new TH1F(name,name,200,0,200);
 334     fOutput->Add(fHistJetPtBias[icent]);
 335
 336     snprintf(name,200,"fHistJetPtTT_%i",icent);
 337     fHistJetPtTT[icent] = new TH1F(name,name,200,0,200);
 338     fOutput->Add(fHistJetPtTT[icent]);
 339
 340     snprintf(name,200,"fHistAreavsRawPt_%i",icent);
 341     fHistAreavsRawPt[icent] = new TH2F(name,name,100,0,1,200,0,200);
 342     //fOutput->Add(fHistAreavsRawPt[icent]);
 343
 344     for(Int_t iptjet = 0; iptjet<5; ++iptjet){
 345       for(Int_t ieta = 0; ieta<3; ++ieta){
 346         snprintf(name,200,"fHistJetH_%i_%i_%i",icent,iptjet,ieta);
 347         fHistJetH[icent][iptjet][ieta]=new TH2F(name,name,64,-0.5*TMath::Pi(),1.5*TMath::Pi(),300,0,30);
 348         fOutput->Add(fHistJetH[icent][iptjet][ieta]);
 349
 350         snprintf(name,200,"fHistJetHBias_%i_%i_%i",icent,iptjet,ieta);
 351         fHistJetHBias[icent][iptjet][ieta]=new TH2F(name,name,64,-0.5*TMath::Pi(),1.5*TMath::Pi(),300,0,30);
 352         fOutput->Add(fHistJetHBias[icent][iptjet][ieta]);
 353
 354         snprintf(name,200,"fHistJetHTT_%i_%i_%i",icent,iptjet,ieta);
 355         fHistJetHTT[icent][iptjet][ieta]=new TH2F(name,name,64,-0.5*TMath::Pi(),1.5*TMath::Pi(),300,0,30);
 356         fOutput->Add(fHistJetHTT[icent][iptjet][ieta]);
 357       } // end of eta loop
 358     } // end of pt-jet loop
 359   } // end of centrality loop
 360
 361   // strings for titles
 362   TString name1;
 363   TString title1;
 364
 365   for (Int_t i = 0;i<6;++i){
 366     name1 = TString(Form("JetPtvsTrackPt_%i",i));
 367     title1 = TString(Form("Jet p_{T} vs Leading Track p_{T} cent bin %i",i));
 368     fHistJetPtvsTrackPt[i] = new TH2F(name1,title1,250,-50,200,250,0,50);
 369     fOutput->Add(fHistJetPtvsTrackPt[i]);
 370
 371     name1 = TString(Form("RawJetPtvsTrackPt_%i",i));
 372     title1 = TString(Form("Raw Jet p_{T} vs Leading Track p_{T} cent bin %i",i));
 373     fHistRawJetPtvsTrackPt[i] = new TH2F(name1,title1,250,-50,200,250,0,50);
 374     fOutput->Add(fHistRawJetPtvsTrackPt[i]);
 375
 376     name1 = TString(Form("TrackPt_%i",i));
 377     title1 = TString(Form("Track p_{T} cent bin %i",i));
 378     fHistTrackPt[i] = new TH1F(name1,title1,1000,0,100); // up to 200?
 379     fOutput->Add(fHistTrackPt[i]);
 380
 381     name1 = TString(Form("EP0_%i",i));
 382     title1 = TString(Form("EP VZero cent bin %i",i));
 383     fHistEP0[i] = new TH1F(name1,title1,144,-TMath::Pi(),TMath::Pi());
 384     fOutput->Add(fHistEP0[i]);
 385
 386 /*
 387     name1 = TString(Form("EP0A_%i",i));
 388     title1 = TString(Form("EP VZero cent bin %i",i));
 389     fHistEP0A[i] = new TH1F(name1,title1,144,-TMath::Pi(),TMath::Pi());
 390     fOutput->Add(fHistEP0A[i]);
 391
 392     name1 = TString(Form("EP0C_%i",i));
 393     title1 = TString(Form("EP VZero cent bin %i",i));
 394     fHistEP0C[i] = new TH1F(name1,title1,144,-TMath::Pi(),TMath::Pi());
 395     fOutput->Add(fHistEP0C[i]);
 396
 397     name1 = TString(Form("EPAvsC_%i",i));
 398     title1 = TString(Form("EP VZero cent bin %i",i));
 399     fHistEPAvsC[i] = new TH2F(name1,title1,144,-TMath::Pi(),TMath::Pi(),144,-TMath::Pi(),TMath::Pi());
 400     fOutput->Add(fHistEPAvsC[i]);
 401 */
 402
 403     name1 = TString(Form("JetPtvsdEP_%i",i));
 404     title1 = TString(Form("Jet p_{T} vs #DeltaEP cent bin %i",i));
 405     fHistJetPtvsdEP[i] = new TH2F(name1,title1,250,-50,200,288,-2*TMath::Pi(),2*TMath::Pi());
 406     fOutput->Add(fHistJetPtvsdEP[i]);
 407
 408     name1 = TString(Form("JetPtvsdEPBias_%i",i));
 409     title1 = TString(Form("Bias Jet p_{T} vs #DeltaEP cent bin %i",i));
 410     fHistJetPtvsdEPBias[i] = new TH2F(name1,title1,250,-50,200,288,-2*TMath::Pi(),2*TMath::Pi());
 411     fOutput->Add(fHistJetPtvsdEPBias[i]);
 412
 413     name1 = TString(Form("RhovsdEP_%i",i));
 414     title1 = TString(Form("#rho vs #DeltaEP cent bin %i",i));
 415     fHistRhovsdEP[i] = new TH2F(name1,title1,500,0,500,288,-2*TMath::Pi(),2*TMath::Pi());
 416     fOutput->Add(fHistRhovsdEP[i]);
 417
 418     name1 = TString(Form("JetEtaPhiPt_%i",i));
 419     title1 = TString(Form("Jet #eta-#phi p_{T} cent bin %i",i));
 420     fHistJetEtaPhiPt[i] = new TH3F(name1,title1,250,-50,200,100,-1,1,64,-3.2,3.2);
 421     fOutput->Add(fHistJetEtaPhiPt[i]);
 422
 423     name1 = TString(Form("JetEtaPhiPtBias_%i",i));
 424     title1 = TString(Form("Jet #eta-#phi p_{T} Bias cent bin %i",i));
 425     fHistJetEtaPhiPtBias[i] = new TH3F(name1,title1,250,-50,200,100,-1,1,64,-3.2,3.2);
 426     fOutput->Add(fHistJetEtaPhiPtBias[i]);
 427
 428     name1 = TString(Form("JetPtArea_%i",i));
 429     title1 = TString(Form("Jet p_{T} Area cent bin %i",i));
 430     fHistJetPtArea[i] = new TH2F(name1,title1,250,-50,200,100,0,1);
 431     fOutput->Add(fHistJetPtArea[i]);
 432
 433     name1 = TString(Form("JetPtAreaBias_%i",i));
 434     title1 = TString(Form("Jet p_{T} Area Bias cent bin %i",i));
 435     fHistJetPtAreaBias[i] = new TH2F(name1,title1,250,-50,200,100,0,1);
 436     fOutput->Add(fHistJetPtAreaBias[i]);
 437
 438     name1 = TString(Form("JetPtNcon_%i",i));
 439     title1 = TString(Form("Jet p_{T} Ncon cent bin %i",i));
 440     fHistJetPtNcon[i] = new TH2F(name1,title1,250,-50,200,100,0,2000);
 441     //fOutput->Add(fHistJetPtNcon[i]);
 442
 443     name1 = TString(Form("JetPtNconBias_%i",i));
 444     title1 = TString(Form("Jet p_{T} NconBias cent bin %i",i));
 445     fHistJetPtNconBias[i] = new TH2F(name1,title1,250,-50,200,100,0,2000);
 446     fOutput->Add(fHistJetPtNconBias[i]);
 447
 448     name1 = TString(Form("JetPtNconCh_%i",i));
 449     title1 = TString(Form("Jet p_{T} NconCh cent bin %i",i));
 450     fHistJetPtNconCh[i] = new TH2F(name1,title1,250,-50,200,100,0,2000);
 451     //fOutput->Add(fHistJetPtNconCh[i]);
 452
 453     name1 = TString(Form("JetPtNconBiasCh_%i",i));
 454     title1 = TString(Form("Jet p_{T} NconBiasCh cent bin %i",i));
 455     fHistJetPtNconBiasCh[i] = new TH2F(name1,title1,250,-50,200,100,0,2000);
 456     //fOutput->Add(fHistJetPtNconBiasCh[i]);
 457
 458     name1 = TString(Form("JetPtNconEm_%i",i));
 459     title1 = TString(Form("Jet p_{T} NconEm cent bin %i",i));
 460     fHistJetPtNconEm[i] = new TH2F(name1,title1,250,-50,200,100,0,2000);
 461     //fOutput->Add(fHistJetPtNconEm[i]);
 462
 463     name1 = TString(Form("JetPtNconBiasEm_%i",i));
 464     title1 = TString(Form("Jet p_{T} NconBiasEm cent bin %i",i));
 465     fHistJetPtNconBiasEm[i] = new TH2F(name1,title1,250,-50,200,100,0,2000);
 466     //fOutput->Add(fHistJetPtNconBiasEm[i]);
 467
 468     name1 = TString(Form("JetHaddPhiINcent_%i",i));
 469     title1 = TString(Form("Jet Hadron #Delta#varphi Distribution IN PLANE cent bin %i",i));
 470     fHistJetHaddPhiINcent[i] = new TH1F(name1,title1,128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
 471     fOutput->Add(fHistJetHaddPhiINcent[i]);
 472
 473     name1 = TString(Form("JetHaddPhiOUTcent_%i",i));
 474     title1 = TString(Form("Jet Hadron #Delta#varphi Distribution OUT PLANE cent bin %i",i));
 475     fHistJetHaddPhiOUTcent[i] = new TH1F(name1,title1,128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
 476     fOutput->Add(fHistJetHaddPhiOUTcent[i]);
 477
 478     name1 = TString(Form("JetHaddPhiMIDcent_%i",i));
 479     title1 = TString(Form("Jet Hadron #Delta#varphi Distribution MIDDLE of PLANE cent bin %i",i));
 480     fHistJetHaddPhiMIDcent[i] = new TH1F(name1,title1,128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
 481     fOutput->Add(fHistJetHaddPhiMIDcent[i]);
 482   }
 483
 484   // add histo's to output list
 485   fOutput->Add(fHistTPCdEdX);
 486   fOutput->Add(fHistITSsignal);
 487 //  fOutput->Add(fHistTOFsignal);
 488   fOutput->Add(fHistRhovsCent);
 489   fOutput->Add(fHistNjetvsCent);
 490   fOutput->Add(fHistTrackPtallcent);
 491   fOutput->Add(fHistJetEtaPhi);
 492   //fOutput->Add(fHistTrackEtaPhi);
 493   fOutput->Add(fHistJetHEtaPhi);
 494   fOutput->Add(fHistCentrality);
 495   fOutput->Add(fHistZvtx);
 496   fOutput->Add(fHistJetPhi);
 497   fOutput->Add(fHistTrackPhi);
 498   fOutput->Add(fHistJetHaddPhiIN);
 499   fOutput->Add(fHistJetHaddPhiOUT);
 500   fOutput->Add(fHistJetHaddPhiMID);
 501   fOutput->Add(fHistJetHaddPhiBias);
 502   fOutput->Add(fHistJetHaddPhiINBias);
 503   fOutput->Add(fHistJetHaddPhiOUTBias);
 504   fOutput->Add(fHistJetHaddPhiMIDBias);
 505   fOutput->Add(fHistSEphieta);
 506   fOutput->Add(fHistMEphieta);
 507   fOutput->Add(fHistJetHaddPHI);
 508
 509   // set up jet-hadron sparse
 510   UInt_t bitcoded = 0; // bit coded, see GetDimParams() below
 511   UInt_t cifras = 0;
 512   UInt_t bitcode = 0;  // bit coded, see GetDimParamsPID() below
 513   bitcoded = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<7 | 1<<8 | 1<<9;
 514   fhnJH = NewTHnSparseF("fhnJH", bitcoded);
 515   fhnJH->Sumw2();
 516 ///  fOutput->Add(fhnJH);
 517
 518   // ****************************** PID *****************************************************
 519   // set up PID handler
 520   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 521   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 522   if(!inputHandler) {
 523         AliError("Input handler needed");
 524         return;
 525   }
 526
 527   // PID response object
 528   fPIDResponse = inputHandler->GetPIDResponse();
 529   if (!fPIDResponse) {
 530         AliError("PIDResponse object was not created");
 531         return;
 532   }
 533
 534   // ****************************************************************************************
 535 /*
 536   // for pp we need mult bins for event mixing. Create binning here, to also make a histogram from it
 537   Int_t nCentralityBins  = 8;
 538   Double_t centralityBins[9] = {0.0, 4., 9, 15, 25, 35, 55, 100.0,500.0};
 539   Double_t centralityBins[nCentralityBins+1];
 540   for(Int_t ic=0; ic<nCentralityBins+1; ic++){
 541     if(ic==nCentralityBins) centralityBins[ic]=500;
 542     else centralityBins[ic]=10.0*ic;
 543   }
 544 */
 545   // setup for Pb-Pb collisions
 546   Int_t nCentralityBins  = 100;
 547   Double_t centralityBins[nCentralityBins];
 548   for(Int_t ic=0; ic<nCentralityBins; ic++){
 549     centralityBins[ic]=1.0*ic;
 550   }
 551
 552   fHistMult = new TH1F("fHistMult","multiplicity",nCentralityBins,centralityBins);
 553   fOutput->Add(fHistMult);
 554
 555   // Event Mixing
 556   Int_t trackDepth = fMixingTracks;
 557   Int_t poolsize   = 1000;  // Maximum number of events, ignored in the present implemented of AliEventPoolManager
 558   Int_t nZvtxBins  = 5+1+5;
 559   Double_t vertexBins[] = { -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10};
 560   Double_t* zvtxbin = vertexBins;
 561   fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centralityBins, nZvtxBins, zvtxbin);
 562
 563   // set up event mixing sparse
 564   if(fDoEventMixing){
 565      cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<7;
 566      fhnMixedEvents = NewTHnSparseF("fhnMixedEvents", cifras);
 567      fhnMixedEvents->Sumw2();
 568      fOutput->Add(fhnMixedEvents);
 569   } // end of do-eventmixing
 570
 571   Double_t xlowjetPT[] = {-50,-45,-40,-35,-30,-25,-20,-18,-16,-14,-12,-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10,12,14,16,18,20,25,30,35,40,45,50,60,70,80,90,100,120,140,160,180,200,250,300,350,400};
 572   Double_t xlowtrPT[] = {0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.25,2.50,2.75,3.0,3.25,3.5,3.75,4.0,4.25,4.50,4.75,5.0,5.5,6.0,6.5,7.0,7.5,8.0,8.5,9.0,9.5,10.0,11.0,12.0,13.0,14.0,15.0,16.0,17.0,18.0,19.0,20.0,22.0,24.0,26.0,28.0,30.0,35.0,40.0,45.0,50.0,60.0,70.0,80.0,90.0,100.0};
 573
 574   // number of bins you tell histogram should be (# in array - 1) because the last bin
 575   // is the right-most edge of the histogram
 576   // i.e. this is for PT and there are 57 numbers (bins) thus we have 56 bins in our histo
 577   Int_t nbinsjetPT = sizeof(xlowjetPT)/sizeof(Double_t) - 1;
 578   Int_t nbinstrPT = sizeof(xlowtrPT)/sizeof(Double_t) - 1;
 579
 580   // set up PID sparse
 581   if(doPID){
 582     bitcode = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<7 | 1<<8 | 1<<9 |
 583               1<<10 | 1<<11 | 1<<12 | 1<<13 | 1<<14 | 1<<15 | 1<<16 | 1<<17 | 1<<18 | 1<<19 |
 584               1<<20 | 1<<21 | 1<<22;
 585     fhnPID = NewTHnSparseFPID("fhnPID", bitcode);
 586     fhnPID->Sumw2();
 587     fOutput->Add(fhnPID);
 588   } // end of do-PID
 589
 590   if(dovarbinTHnSparse){
 591     fhnPID->GetAxis(1)->Set(nbinsjetPT, xlowjetPT);
 592     fhnPID->GetAxis(2)->Set(nbinstrPT, xlowtrPT);
 593   }
 594
 595   // =========== Switch on Sumw2 for all histos ===========
 596   for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
 597     TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
 598     if (h1){
 599       h1->Sumw2();
 600       continue;
 601     }
 602     TH2 *h2 = dynamic_cast<TH2*>(fOutput->At(i));
 603     if (h2){
 604       h2->Sumw2();
 605       continue;
 606     }
 607     TH3 *h3 = dynamic_cast<TH3*>(fOutput->At(i));
 608     if (h3){
 609       h3->Sumw2();
 610       continue;
 611     }
 612     THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
 613     if(hn)hn->Sumw2();
 614   }
 615
 616   PostData(1, fOutput);
 617 }
 618
 619 //_________________________________________________________________________
 620 void AliAnalysisTaskEmcalJetHadEPpid::ExecOnce()
 621 {
 622   AliAnalysisTaskEmcalJet::ExecOnce();
 623 }
 624
 625 //_________________________________________________________________________
 626 Bool_t AliAnalysisTaskEmcalJetHadEPpid::Run()
 627 { // Main loop called for each event
 628   // TEST TEST TEST TEST for OBJECTS!
 629   if(!fLocalRho){
 630     AliError(Form("Couldn't get fLocalRho object, try to get it from Event based on name\n"));
 631     fLocalRho = GetLocalRhoFromEvent(fLocalRhoName);
 632   }
 633   if(!fTracks){
 634     AliError(Form("No fTracks object!!\n"));
 635     return kTRUE;
 636   }
 637   if(!fJets){
 638     AliError(Form("No fJets object!!\n"));
 639     return kTRUE;
 640   }
 641
 642   // what kind of event do we have: AOD or ESD?
 643   if (dynamic_cast<AliAODEvent*>(InputEvent())) useAOD = kTRUE;
 644   else useAOD = kFALSE;
 645
 646   // if we have ESD event, set up ESD object
 647   if(!useAOD){
 648     fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 649     if (!fESD) {
 650       AliError(Form("ERROR: fESD not available\n"));
 651       return kTRUE;
 652     }
 653   }
 654
 655   // if we have AOD event, set up AOD object
 656   if(useAOD){
 657     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 658     if(!fAOD) {
 659       AliError(Form("ERROR: fAOD not available\n"));
 660       return kTRUE;
 661     }
 662   }
 663
 664   // get centrality
 665   Int_t centbin = GetCentBin(fCent);
 666   fHistCentrality->Fill(fCent); // won't be filled in pp collision (Keep this in mind!)
 667
 668   // for pp analyses we will just use the first centrality bin
 669   if (centbin == -1) centbin = 0;
 670
 671   // get vertex information
 672   Double_t fvertex[3]={0,0,0};
 673   InputEvent()->GetPrimaryVertex()->GetXYZ(fvertex);
 674   Double_t zVtx=fvertex[2];
 675   fHistZvtx->Fill(zVtx);
 676
 677   // get z-vertex bin
 678   //Int_t zVbin = GetzVertexBin(zVtx);
 679
 680   // apply zVtx cut
 681   if(fabs(zVtx)>10.0)
 682     return kTRUE;
 683
 684   // apply cut to event on Centrality > 90%
 685   if(fCent>90) return kTRUE;
 686
 687   // create pointer to list of input event
 688   TList *list = InputEvent()->GetList();
 689   if(!list) {
 690     AliError(Form("ERROR: list not attached\n"));
 691     return kFALSE;
 692   }
 693
 694   // initialize TClonesArray pointers to jets and tracks
 695   TClonesArray *jets = 0;
 696   TClonesArray *tracks = 0;
 697
 698   // get Tracks object
 699   tracks = dynamic_cast<TClonesArray*>(list->FindObject(fTracks));
 700   if (!tracks) {
 701     AliError(Form("Pointer to tracks %s == 0", fTracks->GetName()));
 702     return kTRUE;
 703   } // verify existence of tracks
 704
 705   // get Jets object
 706   jets = dynamic_cast<TClonesArray*>(list->FindObject(fJets));
 707   if(!jets){
 708     AliError(Form("Pointer to jets %s == 0", fJets->GetName()));
 709     return kTRUE;
 710   } // verify existence of jets
 711
 712   // get number of jets and tracks
 713   const Int_t Njets = jets->GetEntries();
 714   const Int_t Ntracks = tracks->GetEntries();
 715   if(Ntracks<1)   return kTRUE;
 716   if(Njets<1)     return kTRUE;
 717
 718   fHistMult->Fill(Ntracks);  // fill multiplicity distribution
 719
 720   // initialize track parameters
 721   Int_t iTT=-1;
 722   Double_t ptmax=-10;
 723
 724   // loop over tracks - to get hardest track (highest pt)
 725   for (Int_t iTracks = 0; iTracks < Ntracks; iTracks++){
 726     AliVTrack* track = static_cast<AliVTrack*>(tracks->At(iTracks));
 727     if (!track) {
 728       AliError(Form("Couldn't get VTrack track %d\n", iTracks));
 729       continue;
 730     } // verify existence of tracks
 731
 732     if(TMath::Abs(track->Eta())>0.9) continue;
 733     if(track->Pt()<0.15) continue;
 734     //iCount++;
 735     if(track->Pt()>ptmax){
 736       ptmax=track->Pt();             // max pT track
 737       iTT=iTracks;                   // trigger tracks
 738     } // check if Pt>maxpt
 739
 740     fHistTrackPhi->Fill(track->Phi());
 741     fHistTrackPtallcent->Fill(track->Pt());
 742   } // end of loop over tracks
 743
 744   // find tracks for event (alternative), and loop over
 745   for (int i = 0;i<Ntracks;i++){
 746     AliVParticle *track = static_cast<AliVParticle*>(fTracks->At(i));
 747     if (!track) continue;
 748
 749     if(TMath::Abs(track->Eta())>0.9) continue;
 750     if(track->Pt()<0.15) continue;
 751
 752     fHistTrackPt[centbin]->Fill(track->Pt());
 753   }
 754
 755   // event plane histograms filled
 756   fHistEP0[centbin]->Fill(fEPV0);
 757   //fHistEP0A[centbin]->Fill(fEPV0A);
 758   //fHistEP0C[centbin]->Fill(fEPV0C);
 759   //fHistEPAvsC[centbin]->Fill(fEPV0A,fEPV0C);
 760
 761   // get rho from event and fill relative histo's
 762   fRho = GetRhoFromEvent(fRhoName);
 763   fRhoVal = fRho->GetVal();
 764   fHistRhovsCent->Fill(fCent,fRhoVal);        // Global Rho vs Centrality
 765   fHistRhovsdEP[centbin]->Fill(fRhoVal,fEPV0); // Global Rho vs delta Event Plane angle
 766
 767   // initialize jet parameters
 768   Int_t ijethi=-1;
 769   Double_t highestjetpt=0.0;
 770   Int_t passedTTcut=0;
 771   Int_t NjetAcc = 0;
 772
 773   // loop over jets in an event - to find highest jet pT and apply some cuts
 774   for (Int_t ijet = 0; ijet < Njets; ijet++){
 775     // get our jets
 776     AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(ijet));
 777     if (!jet) continue;
 778
 779     // cuts on jets
 780     if ((jet->Phi()<fPhimin)||(jet->Phi()>fPhimax))
 781       continue;
 782     if ((jet->Eta()<fEtamin)||(jet->Eta()>fEtamax))
 783       continue;
 784     fHistAreavsRawPt[centbin]->Fill(jet->Pt(),jet->Area());
 785     if (jet->Area()<fAreacut)
 786       continue;
 787     //prevents 0 area jets from sneaking by when area cut == 0
 788     if (jet->Area()==0)
 789       continue;
 790     //exclude jets with extremely high pt tracks which are likely misreconstructed
 791     if(jet->MaxTrackPt()>100)
 792       continue;
 793     if(!AcceptJet(jet))
 794       continue;
 795
 796     NjetAcc++;                     // # of accepted jets
 797     fHistJetPhi->Fill(jet->Phi()); // Jet Phi histogram (filled)
 798
 799     // get highest Pt jet in event
 800     if(highestjetpt<jet->Pt()){
 801       ijethi=ijet;
 802       highestjetpt=jet->Pt();
 803     }
 804   } // end of looping over jets
 805
 806   // accepted jets
 807   fHistNjetvsCent->Fill(fCent,NjetAcc);
 808   Int_t NJETAcc = 0;
 809
 810   // loop over jets in event and make appropriate cuts
 811   for (Int_t ijet = 0; ijet < Njets; ++ijet) {
 812      AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(ijet));
 813      if (!jet)  // see if we have a jet
 814        continue;
 815      if(!AcceptJet(jet))             // do we accept this jet?
 816        continue;
 817
 818 //     Double_t leadjet=0;
 819 //     if (ijet==ijethi) leadjet=1;
 820
 821      // cuts on jet
 822      if ((jet->Phi()<fPhimin)||(jet->Phi()>fPhimax))
 823         continue;
 824      if ((jet->Eta()<fEtamin)||(jet->Eta()>fEtamax))
 825         continue;
 826      if (jet->Area()==0) // make sure jet has an area
 827        continue;
 828      if (jet->Pt()<0.1) // (should probably be higher..., but makes a cut on jet pT)
 829        continue;
 830      if (jet->MaxTrackPt()>100) // elminates fake tracks
 831        continue;
 832      if (!AcceptJet(jet)) // sees if jet is accepted
 833        continue;
 834      //  jets.push_back(jet);
 835
 836      // check on lead jet
 837      Double_t leadjet=0;
 838      if (ijet==ijethi) leadjet=1;
 839
 840      // initialize and calculate various parameters: pt, eta, phi, rho, etc...
 841      Double_t jetphi = jet->Phi();      // phi of jet
 842      NJETAcc++;   // # accepted jets
 843      fLocalRhoVal = fLocalRho->GetLocalVal(jetphi, GetJetRadius(0)); //fmyJetRadius);   // get local rho value
 844      Double_t jeteta = jet->Eta();     // ETA of jet
 845 //     Double_t jetptraw = jet->Pt();    // raw pT of jet
 846      Double_t jetPt = -500;            // initialize corr jet pt
 847 //     jetPt = jet->Pt();
 848 //     jetPt = jet->Pt()-jet->Area()*fRhoVal;  // corrected pT of jet from rho value
 849      jetPt = jet->Pt()-jet->Area()*fLocalRhoVal; // corrected pT of jet using Rho modulated for V2 and V3
 850      Float_t dEP = -500;                    // initialize angle between jet and event plane
 851      //Int_t ieta = -500;             // initialize jet eta bin
 852      //ieta = GetEtaBin(jeteta);         // bin of eta
 853      dEP = RelativeEPJET(jetphi,fEPV0); // angle betweeen jet and event plane
 854
 855      // make histo's
 856      fHistJetPtvsTrackPt[centbin]->Fill(jetPt,jet->MaxTrackPt());
 857      fHistRawJetPtvsTrackPt[centbin]->Fill(jet->Pt(),jet->MaxTrackPt());
 858 //     fHistJetPtvsdEP[centbin]->Fill(jetPt,RelativePhi((fEPV0+TMath::Pi()),jet->Phi()));
 859      fHistJetPtvsdEP[centbin]->Fill(jetPt, dEP);
 860      fHistJetEtaPhiPt[centbin]->Fill(jetPt,jet->Eta(),jet->Phi());
 861      fHistJetPtArea[centbin]->Fill(jetPt,jet->Area());
 862      fHistJetPtNcon[centbin]->Fill(jetPt,jet->GetNumberOfConstituents());
 863      fHistJetPtNconCh[centbin]->Fill(jetPt,jet->GetNumberOfTracks());
 864      fHistJetPtNconEm[centbin]->Fill(jetPt,jet->GetNumberOfClusters());
 865
 866      // make histo's with BIAS applied
 867      if (jet->MaxTrackPt()>fTrkBias){
 868        fHistJetPtvsdEPBias[centbin]->Fill(jetPt, dEP);
 869        fHistJetEtaPhiPtBias[centbin]->Fill(jetPt,jet->Eta(),jet->Phi());
 870        fHistJetPtAreaBias[centbin]->Fill(jetPt,jet->Area());
 871        fHistJetPtNconBias[centbin]->Fill(jetPt,jet->GetNumberOfConstituents());
 872        fHistJetPtNconBiasCh[centbin]->Fill(jetPt,jet->GetNumberOfTracks());
 873        fHistJetPtNconBiasEm[centbin]->Fill(jetPt,jet->GetNumberOfClusters());
 874      }
 875
 876     // fill some histo's
 877     fHistJetPt[centbin]->Fill(jet->Pt());  // fill #jets vs pT histo
 878     fHistJetEtaPhi->Fill(jet->Eta(),jet->Phi());     // fill jet eta-phi distribution histo
 879     //fHistDeltaPtvsArea->Fill(jetPt,jet->Area());
 880
 881     if(leadjet && centbin==0) fHistJetPt[centbin+1]->Fill(jet->Pt());
 882     if ((jet->MaxTrackPt()>fTrkBias) || (jet->MaxClusterPt()>fClusBias)){
 883       fHistJetPtBias[centbin]->Fill(jet->Pt());
 884       if(leadjet && centbin==0) fHistJetPtBias[centbin+1]->Fill(jet->Pt());
 885     }  // end of MaxTrackPt>ftrkBias or maxclusterPt>fclusBias
 886
 887     // do we have trigger tracks
 888     if(iTT>0){
 889       AliVTrack* TT = static_cast<AliVTrack*>(tracks->At(iTT));
 890       if(TMath::Abs(jet->Phi()-TT->Phi()-TMath::Pi())<0.6) passedTTcut=1;
 891       else passedTTcut=0;
 892     } // end of check on iTT > 0
 893
 894     if(passedTTcut)  fHistJetPtTT[centbin]->Fill(jet->Pt());
 895
 896     // cut on HIGHEST jet pt in event (15 GeV default)
 897     if (highestjetpt>fJetPtcut) {
 898
 899       // loop over all track for an event containing a jet with a pT>fJetPtCut  (15)GeV
 900       for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
 901         //AliVTrack* track = static_cast<AliVTrack*>(tracks->At(iTracks));
 902         // check on even type to get object for PID info
 903         //char eventtype[10];
 904
 905         AliVTrack *track = fAOD->GetTrack(iTracks);
 906         if(!track) track = fESD->GetTrack(iTracks);
 907 //      AliAODTrack* track = fAOD->GetTrack(iTracks);
 908 //         //sprintf(eventtype,"fAOD");
 909 //         if(!track) AliESDtrack* track = fESD->GetTrack(iTracks);
 910 //         //sprintf(eventtype,"fESD");
 911
 912         if (!track) {
 913           AliError(Form("Couldn't get AliESDtrack %d\n", iTracks));
 914           continue;
 915         }
 916
 917 //     AliESDtrack *track = GetAcceptTrack(esdtrack);
 918 //     if(!track) continue;
 919
 920         // apply track cuts
 921         if(TMath::Abs(track->Eta())>fTrkEta)
 922           continue;
 923         //fHistTrackPt->Fill(track->Pt());  // fill track pT histogram
 924         if (track->Pt()<0.15)             // make sure track pT > 0.15 GeV
 925           continue;
 926
 927         // calculate and get some track parameters
 928         track->GetTPCsignal();            // get the TPC signal of the track
 929         Double_t tracketa=track->Eta();   // eta of track
 930         Double_t deta=tracketa-jeteta;    // dETA between track and jet
 931         Int_t ieta=GetEtaBin(deta);       // bin of eta
 932         if(ieta<0) continue;
 933
 934         // dPHI between jet and hadron
 935         Double_t dphijh = RelativePhi(jet->Phi(), track->Phi()); // angle between jet and hadron
 936
 937         // jet pt bins
 938         //Int_t iptjet=-1;                  // initialize jet pT bin
 939 //        iptjet=GetpTjetBin(jetPt);        // bin of jet pT
 940         Int_t iptjet=GetpTjetBin(jet->Pt());
 941         if(iptjet<0) continue;
 942         Double_t dR=sqrt(deta*deta+dphijh*dphijh);                   // difference of R between jet and hadron track
 943
 944         // fill some jet-hadron histo's
 945         fHistJetH[centbin][iptjet][ieta]->Fill(dphijh,track->Pt());  // fill jet-hadron dPHI--track pT distribution
 946         fHistJetHEtaPhi->Fill(deta,dphijh);                          // fill jet-hadron  eta--phi distribution
 947         fHistJetHaddPHI->Fill(dphijh);
 948
 949         // does our max track or cluster pass the bias?
 950         if ((jet->MaxTrackPt()>fTrkBias) || (jet->MaxClusterPt()>fClusBias)){
 951           // set up and fill Jet-Hadron THnSparse
 952           Double_t triggerEntries[9] = {fCent,jetPt,track->Pt(),dR,deta,dphijh,0,leadjet,zVtx};
 953           fhnJH->Fill(triggerEntries);    // fill Sparse Histo with trigger entries
 954
 955           // fill histo's
 956           fHistSEphieta->Fill(dphijh, deta);
 957           fHistJetHBias[centbin][iptjet][ieta]->Fill(dphijh,track->Pt());
 958           fHistJetHaddPhiBias->Fill(dphijh);
 959
 960           // in plane and out of plane histo's
 961           if( dEP>0 && dEP<=(TMath::Pi()/6) ){
 962              // we are IN plane
 963              fHistJetHaddPhiINBias->Fill(dphijh);
 964           }else if( dEP>(TMath::Pi()/3) && dEP<=(TMath::Pi()/2) ){
 965              // we are OUT of PLANE
 966              fHistJetHaddPhiOUTBias->Fill(dphijh);
 967           }else if( dEP>(TMath::Pi()/6) && dEP<=(TMath::Pi()/3) ){
 968              // we are in middle of plane
 969              fHistJetHaddPhiMIDBias->Fill(dphijh);
 970           }
 971         } // end of check maxtrackpt>ftrackbias or maxclusterpt>fclustbias
 972
 973         if(passedTTcut) fHistJetHTT[centbin][iptjet][ieta]->Fill(dphijh,track->Pt());
 974
 975 // **************************************************************************************************************
 976 // *********************************** PID **********************************************************************
 977 // **************************************************************************************************************
 978     if(doPID){
 979       // data information for pid
 980       Double_t eta, pt, EovP, fClsE = -1.;
 981
 982       EovP = 0;
 983       fClsE = 0;
 984
 985       // nSigma of particles in TPC, TOF, and ITS
 986       Double_t nSigmaPion_TPC, nSigmaProton_TPC, nSigmaKaon_TPC = -1.;
 987       Double_t nSigmaPion_TOF, nSigmaProton_TOF, nSigmaKaon_TOF = -1.;
 988       Double_t nSigmaPion_ITS, nSigmaProton_ITS, nSigmaKaon_ITS = -1.;
 989
 990       // misc quantities
 991       Double_t dEdx, ITSsig, TOFsig  = -1.;
 992       Double_t charge     = -99;
 993       //Double_t mom=track->GetTPCmomentum();
 994
 995       // get parameters of track
 996       charge = track->Charge();    // charge of track
 997       eta    = track->Eta();       // ETA of track
 998       pt     = track->Pt();        // pT of track
 999
1000       // extra attempt
1001       AliVEvent *vevent=InputEvent();
1002       if (!vevent||!fPIDResponse) return kTRUE; // just return, maybe put at beginning
1003
1004       // get impact parameters and detector signals
1005       dEdx = track->GetTPCsignal();
1006       ITSsig = track->GetITSsignal();
1007       TOFsig = track->GetTOFsignal();
1008
1009       // fill detector signal histograms
1010       fHistTPCdEdX->Fill(pt, dEdx);
1011       fHistITSsignal->Fill(pt, ITSsig);
1012 //     fHistTOFsignal->Fill(pt, TOFsig);
1013
1014       // TPC nSigma's
1015       nSigmaPion_TPC = fPIDResponse->NumberOfSigmasTPC(track,AliPID::kPion);
1016       nSigmaKaon_TPC = fPIDResponse->NumberOfSigmasTPC(track,AliPID::kKaon);
1017       nSigmaProton_TPC = fPIDResponse->NumberOfSigmasTPC(track,AliPID::kProton);
1018
1019       // TOF nSigma's
1020       nSigmaPion_TOF = fPIDResponse->NumberOfSigmasTOF(track,AliPID::kPion);
1021       nSigmaKaon_TOF = fPIDResponse->NumberOfSigmasTOF(track,AliPID::kKaon);
1022       nSigmaProton_TOF = fPIDResponse->NumberOfSigmasTOF(track,AliPID::kProton);
1023
1024       // ITS nSigma's
1025       nSigmaPion_ITS = fPIDResponse->NumberOfSigmasITS(track,AliPID::kPion);
1026       nSigmaKaon_ITS = fPIDResponse->NumberOfSigmasITS(track,AliPID::kKaon);
1027       nSigmaProton_ITS = fPIDResponse->NumberOfSigmasITS(track,AliPID::kProton);
1028
1029       // Tests to PID pions, kaons, and protons,          (default is undentified tracks)
1030       //Double_t nPIDtpc, nPIDits, nPIDtof = 0;
1031       Double_t  nPID = -10;
1032
1033       // check track has pT < 0.900 GeV  - use TPC pid
1034       if (pt<0.900 && dEdx>0) {
1035         nPIDtpc = 4;
1036
1037         // PION check - TPC
1038         if (TMath::Abs(nSigmaPion_TPC)<2 && TMath::Abs(nSigmaKaon_TPC)>2 && TMath::Abs(nSigmaProton_TPC)>2 ){
1039           isPItpc = kTRUE;
1040           nPIDtpc = 1;
1041           nPID=0.5;
1042         }else isPItpc = kFALSE;
1043
1044         // KAON check - TPC
1045         if (TMath::Abs(nSigmaKaon_TPC)<2 && TMath::Abs(nSigmaPion_TPC)>3 && TMath::Abs(nSigmaProton_TPC)>2 ){
1046           isKtpc = kTRUE;
1047           nPIDtpc = 2;
1048           nPID=1.5;
1049         }else isKtpc = kFALSE;
1050
1051         // PROTON check - TPC
1052         if (TMath::Abs(nSigmaProton_TPC)<2 && TMath::Abs(nSigmaPion_TPC)>3 && TMath::Abs(nSigmaKaon_TPC)>2 ){
1053           isPtpc = kTRUE;
1054           nPIDtpc = 3;
1055           nPID=2.5;
1056         }else isPtpc = kFALSE;
1057       }  // cut on track pT for TPC
1058
1059       // check track has pT < 0.500 GeV - use ITS pid
1060       if (pt<0.500 && ITSsig>0) {
1061         nPIDits = 4;
1062
1063         // PION check - ITS
1064         if (TMath::Abs(nSigmaPion_ITS)<2 && TMath::Abs(nSigmaKaon_ITS)>2 && TMath::Abs(nSigmaProton_ITS)>2 ){
1065           isPIits = kTRUE;
1066           nPIDits = 1;
1067           nPID=4.5;
1068         }else isPIits = kFALSE;
1069
1070         // KAON check - ITS
1071         if (TMath::Abs(nSigmaKaon_ITS)<2 && TMath::Abs(nSigmaPion_ITS)>3 && TMath::Abs(nSigmaProton_ITS)>2 ){
1072           isKits = kTRUE;
1073           nPIDits = 2;
1074           nPID=5.5;
1075         }else isKits = kFALSE;
1076
1077         // PROTON check - ITS
1078         if (TMath::Abs(nSigmaProton_ITS)<2 && TMath::Abs(nSigmaPion_ITS)>3 && TMath::Abs(nSigmaKaon_ITS)>2 ){
1079           isPits = kTRUE;
1080           nPIDits = 3;
1081           nPID=6.5;
1082         }else isPits = kFALSE;
1083       }  // cut on track pT for ITS
1084
1085       // check track has 0.900 GeV < pT < 2.500 GeV - use TOF pid
1086       if (pt>0.900 && pt<2.500 && TOFsig>0) {
1087         nPIDtof = 4;
1088
1089         // PION check - TOF
1090         if (TMath::Abs(nSigmaPion_TOF)<2 && TMath::Abs(nSigmaKaon_TOF)>2 && TMath::Abs(nSigmaProton_TOF)>2 ){
1091           isPItof = kTRUE;
1092           nPIDtof = 1;
1093           nPID=8.5;
1094         }else isPItof = kFALSE;
1095
1096         // KAON check - TOF
1097         if (TMath::Abs(nSigmaKaon_TOF)<2 && TMath::Abs(nSigmaPion_TOF)>3 && TMath::Abs(nSigmaProton_TOF)>2 ){
1098           isKtof = kTRUE;
1099           nPIDtof = 2;
1100           nPID=9.5;
1101         }else isKtof = kFALSE;
1102
1103         // PROTON check - TOF
1104         if (TMath::Abs(nSigmaProton_TOF)<2 && TMath::Abs(nSigmaPion_TOF)>3 && TMath::Abs(nSigmaKaon_TOF)>2 ){
1105           isPtof = kTRUE;
1106           nPIDtof = 3;
1107           nPID=10.5;
1108         }else isPtof = kFALSE;
1109       }  // cut on track pT for TOF
1110
1111          if (nPIDtpc == 0) nPID = -0.5;
1112          if (nPIDits == 0) nPID = 3.5;
1113          if (nPIDtof == 0) nPID = 7.5;
1114          fHistPID->Fill(nPID);
1115
1116           // PID sparse getting filled
1117           Double_t nontrig_tracks_Entries[23] = {fCent,jetPt,pt,charge,eta,deta,dphijh,leadjet,zVtx,EovP,fClsE,
1118                                       nSigmaPion_TPC, nSigmaProton_TPC, nSigmaKaon_TPC,  //nSig in TPC
1119                                       nSigmaPion_ITS, nSigmaProton_ITS, nSigmaKaon_ITS,  //nSig in ITS
1120                                       nSigmaPion_TOF, nSigmaProton_TOF, nSigmaKaon_TOF,  //nSig in TOF
1121                                       nPIDtpc, nPIDits, nPIDtof                          //PID label for each detector
1122                                       }; //array for PID sparse
1123           fhnPID->Fill(nontrig_tracks_Entries);   // fill Sparse histo of PID tracks
1124         } // end of doPID check
1125
1126         // get track pt bin
1127         Int_t itrackpt = -500;              // initialize track pT bin
1128         itrackpt = GetpTtrackBin(track->Pt());
1129
1130         // all tracks: jet hadron correlations in hadron pt bins
1131         fHistJetHadbindPhi[itrackpt]->Fill(dphijh);
1132
1133         // in plane and out of plane jet-hadron histo's
1134         if( dEP>0 && dEP<=(TMath::Pi()/6) ){
1135            // we are IN plane
1136            fHistJetHaddPhiINcent[centbin]->Fill(dphijh);
1137            fHistJetHaddPhiIN->Fill(dphijh);
1138            fHistJetHadbindPhiIN[itrackpt]->Fill(dphijh);
1139            //fHistJetHaddPhiPtcentbinIN[itrackpt][centbin]->Fill(dphijh);
1140         }else if( dEP>(TMath::Pi()/3) && dEP<=(TMath::Pi()/2) ){
1141            // we are OUT of PLANE
1142            fHistJetHaddPhiOUTcent[centbin]->Fill(dphijh);
1143            fHistJetHaddPhiOUT->Fill(dphijh);
1144            fHistJetHadbindPhiOUT[itrackpt]->Fill(dphijh);
1145            //fHistJetHaddPhiPtcentbinOUT[itrackpt][centbin]->Fill(dphijh);
1146         }else if( dEP>(TMath::Pi()/6) && dEP<=(TMath::Pi()/3) ){
1147            // we are in the middle of plane
1148            fHistJetHaddPhiMIDcent[centbin]->Fill(dphijh);
1149            fHistJetHaddPhiMID->Fill(dphijh);
1150            fHistJetHadbindPhiMID[itrackpt]->Fill(dphijh);
1151         }
1152       } // loop over tracks found in event with highest JET pT > 10.0 GeV (change)
1153     } // jet pt cut
1154   } // jet loop
1155
1156 // ***************************************************************************************************************
1157 // ******************************** Event MIXING *****************************************************************
1158 //  TObjArray* tracksClone = CloneAndReduceTrackList(tracks); // TEST
1159
1160   //Prepare to do event mixing
1161   if(fDoEventMixing>0){
1162     // event mixing
1163
1164     // 1. First get an event pool corresponding in mult (cent) and
1165     //    zvertex to the current event. Once initialized, the pool
1166     //    should contain nMix (reduced) events. This routine does not
1167     //    pre-scan the chain. The first several events of every chain
1168     //    will be skipped until the needed pools are filled to the
1169     //    specified depth. If the pool categories are not too rare, this
1170     //    should not be a problem. If they are rare, you could lose
1171     //    statistics.
1172
1173     // 2. Collect the whole pool's content of tracks into one TObjArray
1174     //    (bgTracks), which is effectively a single background super-event.
1175
1176     // 3. The reduced and bgTracks arrays must both be passed into
1177     //    FillCorrelations(). Also nMix should be passed in, so a weight
1178     //    of 1./nMix can be applied.
1179
1180     // mix jets from triggered events with tracks from MB events
1181     // get the trigger bit
1182     // need to change trigger bits between different runs
1183     UInt_t trigger = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
1184     if (trigger==0)  return kTRUE; // return
1185
1186     //Double_t Ntrks=(Double_t)Ntracks*1.0;
1187     //cout<<"Test.. Ntrks: "<<fPoolMgr->GetEventPool(Ntrks);
1188
1189     AliEventPool* pool = fPoolMgr->GetEventPool(fCent, zVtx); // for PbPb? fcent
1190     //AliEventPool* pool = fPoolMgr->GetEventPool(Ntrks, zVtx); // for pp
1191
1192     if (!pool){
1193       AliFatal(Form("No pool found for centrality = %f, zVtx = %f", fCent, zVtx));
1194       //AliFatal(Form("No pool found for multiplicity = %f, zVtx = %f", Ntrks, zVtx));
1195       return kTRUE;
1196     }
1197
1198     // use only jets from EMCal-triggered events (for lhc11a use AliVEvent::kEMC1)
1199 ///    if (trigger & AliVEvent::kEMC1) {
1200     if (trigger & AliVEvent::kEMCEJE) {  // TEST
1201       //check for a trigger jet
1202       // fmixingtrack/10 ??
1203       if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 /*100*/ || pool->GetCurrentNEvents() >= 5) {
1204         // loop over jets (passing cuts?)
1205         for (Int_t ijet = 0; ijet < Njets; ijet++) {
1206           //Double_t leadjet=0;
1207           //if (ijet==ijethi) leadjet=1;
1208
1209           // get jet object
1210           AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(ijet));
1211           if(!AcceptJet(jet))  continue;
1212
1213           // get jet stats
1214           Double_t jetphi = jet->Phi();  // phi of jet
1215           Double_t jetPt = jet->Pt();    // pT of jet
1216
1217           Int_t nMix = pool->GetCurrentNEvents();  // how many particles in pool to mix
1218
1219           // Fill for biased jet triggers only
1220           if ((jet->MaxTrackPt()>fTrkBias) || (jet->MaxClusterPt()>fClusBias)) {
1221             // Fill mixed-event histos here
1222             for (Int_t jMix=0; jMix<nMix; jMix++) {
1223                 TObjArray* bgTracks = pool->GetEvent(jMix);
1224                 const Int_t Nbgtrks = bgTracks->GetEntries();
1225                 for(Int_t ibg=0; ibg<Nbgtrks; ibg++) {
1226                   AliPicoTrack *part = static_cast<AliPicoTrack*>(bgTracks->At(ibg));
1227                   if(!part) continue;
1228
1229                   //Double_t DPhi = jetphi - part->Phi();
1230                   Double_t DEta = part->Eta()-jet->Eta();                // difference in eta
1231                   Double_t DPhi = RelativePhi(jet->Phi(),part->Phi());   // difference in phi
1232                   Float_t dEP = RelativeEPJET(jetphi,fEPV0);
1233                   Double_t charge = part->Charge();
1234
1235                   //Double_t DR=TMath::Sqrt(DPhi*DPhi+DEta*DEta);      // difference in R
1236                   if(DPhi<-0.5*TMath::Pi()) DPhi+=2.*TMath::Pi();
1237                   if(DPhi>3./2.*TMath::Pi()) DPhi-=2.*TMath::Pi();
1238
1239                   Double_t triggerEntries[8] = {fCent,jetPt,part->Pt(),DEta,DPhi,dEP,zVtx, charge}; //array for ME sparse
1240                   fhnMixedEvents->Fill(triggerEntries,1./nMix);   // fill Sparse histo of mixed events
1241
1242                   // *********************
1243                   //create phi distribution from mixed events
1244                   //fHistMEdPHI->Fill(DPhi);
1245                   fHistMEphieta->Fill(DPhi,DEta, 1./nMix);
1246
1247                 } // end of background track loop
1248              } // end of filling mixed-event histo's
1249           } // end of check for biased jet triggers
1250         } // end of jet loop
1251       } // end of check for triggered jet
1252     } //end EMC triggered loop
1253
1254     // use only tracks from MB events (for lhc11a use AliVEvent::kMB)
1255 ///    if (trigger & AliVEvent::kMB) {
1256     if (trigger & AliVEvent::kAnyINT){ // test
1257       // create a list of reduced objects. This speeds up processing and reduces memory consumption for the event pool
1258       TObjArray* tracksClone = CloneAndReduceTrackList(tracks);
1259
1260       // update pool if jet in event or not
1261       pool->UpdatePool(tracksClone);
1262     } // check on track from MB events
1263   } // end of event mixing
1264
1265   // print some stats on the event
1266   event++;
1267
1268   //cout<<"Event #: "<<event<<endl;
1269   //cout<<"# of jets: "<<Njets<<"      Highest jet pt: "<<highestjetpt<<endl;
1270   //cout<<"# tracks: "<<Ntracks<<"      Highest track pt: "<<ptmax<<endl;
1271   //cout<<" =============================================== "<<endl;
1272
1273   return kTRUE;  // used when the function is of type bool
1274 }  // end of RUN
1275
1276 //________________________________________________________________________
1277 Int_t AliAnalysisTaskEmcalJetHadEPpid::GetCentBin(Double_t cent) const
1278 {  // Get centrality bin.
1279   Int_t centbin = -1;
1280   if (cent>=0 && cent<10)       centbin = 0;
1281   else if (cent>=10 && cent<20) centbin = 1;
1282   else if (cent>=20 && cent<30) centbin = 2;
1283   else if (cent>=30 && cent<40) centbin = 3;
1284   else if (cent>=40 && cent<50) centbin = 4;
1285   else if (cent>=50 && cent<90) centbin = 5;
1286
1287   return centbin;
1288 }
1289
1290 //________________________________________________________________________
1291 Float_t AliAnalysisTaskEmcalJetHadEPpid::RelativePhi(Double_t mphi,Double_t vphi) const
1292 { // function to calculate relative PHI
1293   double dphi = mphi-vphi;
1294
1295   // set dphi to operate on adjusted scale
1296   if(dphi<-0.5*TMath::Pi()) dphi+=2.*TMath::Pi();
1297   if(dphi>3./2.*TMath::Pi()) dphi-=2.*TMath::Pi();
1298
1299   // test
1300   if( dphi < -1.*TMath::Pi()/2 || dphi > 3.*TMath::Pi()/2 )
1301     AliWarning(Form("%s: dPHI not in range [-0.5*Pi, 1.5*Pi]!", GetName()));
1302
1303   return dphi; // dphi in [-0.5Pi, 1.5Pi]
1304 }
1305
1306
1307
1308 //_________________________________________________________________________
1309 Float_t AliAnalysisTaskEmcalJetHadEPpid:: RelativeEPJET(Double_t jetAng, Double_t EPAng) const
1310 { // function to calculate angle between jet and EP in the 1st quadrant (0,Pi/2)
1311   Double_t dphi = (EPAng - jetAng);
1312
1313   // ran into trouble with a few dEP<-Pi so trying this...
1314   if( dphi<-1*TMath::Pi() ){
1315     dphi = dphi + 1*TMath::Pi();
1316   } // this assumes we are doing full jets currently
1317
1318   if( (dphi>0) && (dphi<1*TMath::Pi()/2) ){
1319     // Do nothing! we are in quadrant 1
1320   }else if( (dphi>1*TMath::Pi()/2) && (dphi<1*TMath::Pi()) ){
1321     dphi = 1*TMath::Pi() - dphi;
1322   }else if( (dphi<0) && (dphi>-1*TMath::Pi()/2) ){
1323     dphi = fabs(dphi);
1324   }else if( (dphi<-1*TMath::Pi()/2) && (dphi>-1*TMath::Pi()) ){
1325     dphi = dphi + 1*TMath::Pi();
1326   }
1327
1328   // test
1329   if( dphi < 0 || dphi > TMath::Pi()/2 )
1330     AliWarning(Form("%s: dPHI not in range [0, 0.5*Pi]!", GetName()));
1331
1332   return dphi;   // dphi in [0, Pi/2]
1333 }
1334
1335 //Int_t ieta=GetEtaBin(deta);
1336 //________________________________________________________________________
1337 Int_t AliAnalysisTaskEmcalJetHadEPpid::GetEtaBin(Double_t eta) const
1338 {
1339   // Get eta bin for histos.
1340   Int_t etabin = -1;
1341   if (TMath::Abs(eta)<=0.4)                             etabin = 0;
1342   else if (TMath::Abs(eta)>0.4 && TMath::Abs(eta)<0.8)  etabin = 1;
1343   else if (TMath::Abs(eta)>=0.8)                        etabin = 2;
1344
1345   return etabin;
1346 } // end of get-eta-bin
1347
1348 //________________________________________________________________________
1349 Int_t AliAnalysisTaskEmcalJetHadEPpid::GetpTjetBin(Double_t pt) const
1350 {
1351   // Get jet pt  bin for histos.
1352   Int_t ptbin = -1;
1353   if (pt>=15 && pt<20)          ptbin = 0;
1354   else if (pt>=20 && pt<25)     ptbin = 1;
1355   else if (pt>=25 && pt<40)     ptbin = 2;
1356   else if (pt>=40 && pt<60)     ptbin = 3;
1357   else if (pt>=60)              ptbin = 4;
1358
1359   return ptbin;
1360 } // end of get-jet-pt-bin
1361
1362 //________________________________________________________________________
1363 Int_t AliAnalysisTaskEmcalJetHadEPpid::GetpTtrackBin(Double_t pt) const
1364 {
1365   // May need to update bins for future runs... (testing locally)
1366
1367   // Get track pt bin for histos.
1368   Int_t ptbin = -1;
1369   if (pt < 0.5)                 ptbin = 0;
1370   else if (pt>=0.5 && pt<1.0)   ptbin = 1;
1371   else if (pt>=1.0 && pt<1.5)   ptbin = 2;
1372   else if (pt>=1.5 && pt<2.0)   ptbin = 3;
1373   else if (pt>=2.0 && pt<2.5)   ptbin = 4;
1374   else if (pt>=2.5 && pt<3.0)   ptbin = 5;
1375   else if (pt>=3.0 && pt<4.0)   ptbin = 6;
1376   else if (pt>=4.0 && pt<5.0)   ptbin = 7;
1377   else if (pt>=5.0)             ptbin = 8;
1378
1379   return ptbin;
1380 } // end of get-jet-pt-bin
1381
1382
1383 //___________________________________________________________________________
1384 Int_t AliAnalysisTaskEmcalJetHadEPpid::GetzVertexBin(Double_t zVtx) const
1385 {
1386   // get z-vertex bin for histo.
1387   int zVbin= -1;
1388   if (zVtx>=-10 && zVtx<-8)     zVbin = 0;
1389   else if (zVtx>=-8 && zVtx<-6) zVbin = 1;
1390   else if (zVtx>=-6 && zVtx<-4) zVbin = 2;
1391   else if (zVtx>=-4 && zVtx<-2) zVbin = 3;
1392   else if (zVtx>=-2 && zVtx<0)  zVbin = 4;
1393   else if (zVtx>=0 && zVtx<2)   zVbin = 5;
1394   else if (zVtx>=2 && zVtx<4)   zVbin = 6;
1395   else if (zVtx>=4 && zVtx<6)   zVbin = 7;
1396   else if (zVtx>=6 && zVtx<8)   zVbin = 8;
1397   else if (zVtx>=8 && zVtx<10)  zVbin = 9;
1398   else zVbin = 10;
1399
1400   return zVbin;
1401 } // end of get z-vertex bin
1402
1403 //______________________________________________________________________
1404 THnSparse* AliAnalysisTaskEmcalJetHadEPpid::NewTHnSparseF(const char* name, UInt_t entries)
1405 {
1406    // generate new THnSparseF, axes are defined in GetDimParams()
1407    Int_t count = 0;
1408    UInt_t tmp = entries;
1409    while(tmp!=0){
1410       count++;
1411       tmp = tmp &~ -tmp;  // clear lowest bit
1412    }
1413
1414    TString hnTitle(name);
1415    const Int_t dim = count;
1416    Int_t nbins[dim];
1417    Double_t xmin[dim];
1418    Double_t xmax[dim];
1419
1420    Int_t i=0;
1421    Int_t c=0;
1422    while(c<dim && i<32){
1423       if(entries&(1<<i)){
1424
1425          TString label("");
1426          GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
1427          hnTitle += Form(";%s",label.Data());
1428          c++;
1429       }
1430
1431       i++;
1432    }
1433    hnTitle += ";";
1434
1435    return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
1436 } // end of NewTHnSparseF
1437
1438 void AliAnalysisTaskEmcalJetHadEPpid::GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
1439 {
1440    // stores label and binning of axis for THnSparse
1441    const Double_t pi = TMath::Pi();
1442
1443    switch(iEntry){
1444
1445    case 0:
1446       label = "V0 centrality (%)";
1447       nbins = 10;
1448       xmin = 0.;
1449       xmax = 100.;
1450       break;
1451
1452    case 1:
1453       label = "Corrected Jet p_{T}";
1454       nbins = 250;
1455       xmin = -50.;
1456       xmax = 200.;
1457       break;
1458
1459    case 2:
1460       label = "Track p_{T}";
1461       nbins = 100;
1462       xmin = 0.;
1463       xmax = 100.;
1464       break;
1465
1466     case 3:
1467       label = "Relative Eta";
1468       nbins = 64;
1469       xmin = -1.6;
1470       xmax = 1.6;
1471       break;
1472
1473    case 4:
1474       label = "Relative Phi";
1475       nbins = 144;
1476       xmin = -0.5*pi;
1477       xmax = 1.5*pi;
1478       break;
1479
1480   case 5:
1481       label = "Relative angle of Jet and Reaction Plane";
1482       nbins = 144;
1483       xmin = -0.5*pi;
1484       xmax = 1.5*pi;
1485       break;
1486
1487   case 6:
1488       label = "z-vertex";
1489       nbins = 10;
1490       xmin = -10;
1491       xmax =  10;
1492       break;
1493
1494   case 7:
1495       label = "track charge";
1496       nbins =3;
1497       xmin = -1.5;
1498       xmax = 1.5;
1499       break;
1500
1501   case 8: // need to update
1502       label = "leading track";
1503       nbins = 10;
1504       xmin = 0;
1505       xmax = 50;
1506       break;
1507
1508   case 9:  // need to update
1509       label = "leading jet";
1510       nbins = 3;
1511       xmin = -0.5;
1512       xmax = 2.5;
1513       break;
1514    } // end of switch
1515 } // end of getting dim-params
1516
1517 //_________________________________________________
1518 // From CF event mixing code PhiCorrelations
1519 TObjArray* AliAnalysisTaskEmcalJetHadEPpid::CloneAndReduceTrackList(TObjArray* tracks)
1520 {
1521   // clones a track list by using AliPicoTrack which uses much less memory (used for event mixing)
1522   TObjArray* tracksClone = new TObjArray;
1523   tracksClone->SetOwner(kTRUE);
1524
1525   /*
1526   Int_t nTrax = fESD->GetNumberOfTracks();
1527   //cout << "nTrax " << nTrax <<endl;
1528
1529   for (Int_t i = 0; i < nTrax; ++i)
1530     {
1531       AliESDtrack* esdtrack = fESD->GetTrack(i);
1532       if (!esdtrack)
1533         {
1534           AliError(Form("Couldn't get ESD track %d\n", i));
1535           continue;
1536         }
1537
1538       AliESDtrack *particle = GetAcceptTrack(esdtrack);
1539       if(!particle) continue;
1540
1541       cout << "RM Hybrid track : " << i << "  " << particle->Pt() << endl;
1542   */
1543
1544   //cout << "nEntries " << tracks->GetEntriesFast() <<endl;
1545   for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
1546   {
1547     AliVParticle* particle = (AliVParticle*) tracks->At(i);
1548     if(TMath::Abs(particle->Eta())>fTrkEta) continue;
1549     if(particle->Pt()<0.15)continue;
1550
1551     Double_t trackpt=particle->Pt();   // track pT
1552
1553     Int_t trklabel=-1;
1554     trklabel=particle->GetLabel();
1555     //cout << "TRACK_LABEL: " << particle->GetLabel()<<endl;
1556
1557     Int_t hadbin=-1;
1558     if(trackpt<0.5) hadbin=0;
1559     else if(trackpt<1) hadbin=1;
1560     else if(trackpt<2) hadbin=2;
1561     else if(trackpt<3) hadbin=3;
1562     else if(trackpt<5) hadbin=4;
1563     else if(trackpt<8) hadbin=5;
1564     else if(trackpt<20) hadbin=6;
1565
1566 //feb10 comment out
1567     if(hadbin>-1 && trklabel>-1 && trklabel <3) fHistTrackEtaPhi[trklabel][hadbin]->Fill(particle->Eta(),particle->Phi());
1568     if(hadbin>-1) fHistTrackEtaPhi[3][hadbin]->Fill(particle->Eta(),particle->Phi());
1569
1570     //if(hadbin>-1) fHistTrackEtaPhi[hadbin]->Fill(particle->Eta(),particle->Phi());  // TEST
1571
1572     tracksClone->Add(new AliPicoTrack(particle->Pt(), particle->Eta(), particle->Phi(), particle->Charge(), 0, 0, 0, 0));
1573   } // end of looping through tracks
1574
1575   return tracksClone;
1576 }
1577
1578 //____________________________________________________________________________________________
1579 THnSparse* AliAnalysisTaskEmcalJetHadEPpid::NewTHnSparseFPID(const char* name, UInt_t entries)
1580 {
1581    // generate new THnSparseF PID, axes are defined in GetDimParams()
1582    Int_t count = 0;
1583    UInt_t tmp = entries;
1584    while(tmp!=0){
1585       count++;
1586       tmp = tmp &~ -tmp;  // clear lowest bit
1587    }
1588
1589    TString hnTitle(name);
1590    const Int_t dim = count;
1591    Int_t nbins[dim];
1592    Double_t xmin[dim];
1593    Double_t xmax[dim];
1594
1595    Int_t i=0;
1596    Int_t c=0;
1597    while(c<dim && i<32){
1598       if(entries&(1<<i)){
1599
1600          TString label("");
1601          GetDimParamsPID(i, label, nbins[c], xmin[c], xmax[c]);
1602          hnTitle += Form(";%s",label.Data());
1603          c++;
1604       }
1605
1606       i++;
1607    }
1608    hnTitle += ";";
1609
1610    return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
1611 } // end of NewTHnSparseF PID
1612
1613 //________________________________________________________________________________
1614 void AliAnalysisTaskEmcalJetHadEPpid::GetDimParamsPID(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
1615 {
1616    // stores label and binning of axis for THnSparse
1617    const Double_t pi = TMath::Pi();
1618
1619    switch(iEntry){
1620
1621    case 0:
1622       label = "V0 centrality (%)";
1623       nbins = 10;
1624       xmin = 0.;
1625       xmax = 100.;
1626       break;
1627
1628    case 1:
1629       label = "Corrected Jet p_{T}";
1630       nbins = 250;
1631       xmin = -50.;
1632       xmax = 200.;
1633       break;
1634
1635    case 2:
1636       label = "Track p_{T}";
1637       nbins = 500;
1638       xmin = 0.;
1639       xmax = 50.;
1640       break;
1641
1642    case 3:
1643       label = "Charge of Track";
1644       nbins = 3;
1645       xmin = -1.5;
1646       xmax = 1.5;
1647       break;
1648
1649    case 4:
1650       label = "Track Eta";
1651       nbins = 64;
1652       xmin = -1.2;
1653       xmax = 1.2;
1654       break;
1655
1656    case 5:
1657       label = "Relative Eta of Track and Jet";
1658       nbins = 64;
1659       xmin = -1.6;
1660       xmax = 1.6;
1661       break;
1662
1663    case 6:
1664       label = "Relative Phi of Track and Jet";
1665       nbins = 64;
1666       xmin = -0.5*pi;
1667       xmax = 1.5*pi;
1668       break;
1669
1670    case 7:
1671       label = "leading jet";
1672       nbins = 3;
1673       xmin = -.5;
1674       xmax = 2.5;
1675       break;
1676
1677    case 8:
1678       label = "Z-vertex";
1679       nbins = 10;
1680       xmin = -10.;
1681       xmax = 10.;
1682       break;
1683
1684    case 9: // don't need
1685       label = "E/p of track";
1686       nbins = 100;
1687       xmin = 0.;
1688       xmax = 1.;
1689       break;
1690
1691    case 10: // don't need
1692       label = "Cluster Energy";
1693       nbins = 50; //250;
1694       xmin = 0.;
1695       xmax = 50.; //250.;
1696       break;
1697
1698    case 11:
1699       label = "N-Sigma of pions in TPC";
1700       nbins = 100;
1701       xmin = -5.;
1702       xmax = 5.;
1703       break;
1704
1705    case 12:
1706       label = "N-Sigma of protons in TPC";
1707       nbins = 100;
1708       xmin = -5.;
1709       xmax = 5.;
1710       break;
1711
1712    case 13:
1713       label = "N-Sigma of kaons in TPC";
1714       nbins = 100;
1715       xmin = -5.;
1716       xmax = 5.;
1717       break;
1718
1719    case 14:
1720       label = "N-Sigma of pions in ITS";
1721       nbins = 100;
1722       xmin = -5.;
1723       xmax = 5.;
1724       break;
1725
1726    case 15:
1727       label = "N-Sigma of protons in ITS";
1728       nbins = 100;
1729       xmin = -5.;
1730       xmax = 5.;
1731       break;
1732
1733    case 16:
1734       label = "N-Sigma of kaons in ITS";
1735       nbins = 100;
1736       xmin = -5.;
1737       xmax = 5.;
1738       break;
1739
1740    case 17:
1741       label = "N-Sigma of pions in TOF";
1742       nbins = 100;
1743       xmin = -5.;
1744       xmax = 5.;
1745       break;
1746
1747    case 18:
1748       label = "N-Sigma of protons in TOF";
1749       nbins = 100;
1750       xmin = -5.;
1751       xmax = 5.;
1752       break;
1753
1754    case 19:
1755       label = "N-Sigma of kaons in TOF";
1756       nbins = 100;
1757       xmin = -5.;
1758       xmax = 5.;
1759       break;
1760
1761    case 20:
1762       label = "TPC PID determination";
1763       nbins = 4;
1764       xmin = 0.;
1765       xmax = 4.;
1766       break;
1767
1768    case 21:
1769       label = "ITS PID determination";
1770       nbins = 4;
1771       xmin = 0.;
1772       xmax = 4.;
1773       break;
1774
1775    case 22:
1776       label = "TOF PID determination";
1777       nbins = 4;
1778       xmin = 0.;
1779       xmax = 4.;
1780       break;
1781
1782    } // end of switch
1783 } // end of get dimension parameters PID
1784
1785 void AliAnalysisTaskEmcalJetHadEPpid::Terminate(Option_t *)
1786 {
1787   cout<<"#########################"<<endl;
1788   cout<<"#### DONE RUNNING!!! ####"<<endl;
1789   cout<<"#########################"<<endl;
1790
1791 } // end of terminate
1792
1793
1794 //void AliAnalysisTaskEmcalJetHadEPpid::FillAnalysisSummaryHistogram() const
1795 void AliAnalysisTaskEmcalJetHadEPpid::SetfHistPIDcounterLabels(TH1* h) const
1796 {
1797     // fill the analysis summary histrogram, saves all relevant analysis settigns
1798     h->GetXaxis()->SetBinLabel(1, "TPC: Unidentified"); // -0.5
1799     h->GetXaxis()->SetBinLabel(2, "TPC: Pion"); // 0.5
1800     h->GetXaxis()->SetBinLabel(3, "TPC: Kaon"); // 1.5
1801     h->GetXaxis()->SetBinLabel(4, "TPC: Proton"); // 2.5
1802     h->GetXaxis()->SetBinLabel(5, "ITS: Unidentified"); // 3.5
1803     h->GetXaxis()->SetBinLabel(6, "ITS: Pion"); // 4.5
1804     h->GetXaxis()->SetBinLabel(7, "ITS: Kaon"); // 5.5
1805     h->GetXaxis()->SetBinLabel(8, "ITS: Proton"); // 6.5
1806     h->GetXaxis()->SetBinLabel(9, "TOF: Unidentified"); // 7.5
1807     h->GetXaxis()->SetBinLabel(10, "TOF: Pion"); // 8.5
1808     h->GetXaxis()->SetBinLabel(11, "TOF: Kaon"); // 9.5
1809     h->GetXaxis()->SetBinLabel(12, "TOF: Proton"); // 10.5
1810
1811 }