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