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   useAOD(0), fcutType("EMCAL"), doPID(0), doPIDtrackBIAS(0),
  81   doComments(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   DefineInput(0, TChain::Class());
 167   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   useAOD(0), fcutType("EMCAL"), doPID(0), doPIDtrackBIAS(0),
 181   doComments(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   DefineInput(0, TChain::Class());
 267   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     bitcode = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<7 | 1<<8 | 1<<9 |
 629               1<<10 | 1<<11 | 1<<12 | 1<<13 | 1<<14 | 1<<15 | 1<<16 | 1<<17 | 1<<18 | 1<<19 |
 630               1<<20 | 1<<21 | 1<<22;
 631     fhnPID = NewTHnSparseDPID("fhnPID", bitcode);
 632
 633     if(dovarbinTHnSparse){
 634      fhnPID->GetAxis(1)->Set(nbinsjetPT, xlowjetPT);
 635      fhnPID->GetAxis(2)->Set(nbinstrPT, xlowtrPT);
 636     }
 637
 638     fOutput->Add(fhnPID);
 639   } // end of do-PID
 640
 641   // =========== Switch on Sumw2 for all histos ===========
 642   for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
 643     TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
 644     if (h1){
 645       h1->Sumw2();
 646       continue;
 647     }
 648     TH2 *h2 = dynamic_cast<TH2*>(fOutput->At(i));
 649     if (h2){
 650       h2->Sumw2();
 651       continue;
 652     }
 653     TH3 *h3 = dynamic_cast<TH3*>(fOutput->At(i));
 654     if (h3){
 655       h3->Sumw2();
 656       continue;
 657     }
 658     THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
 659     if(hn)hn->Sumw2();
 660   }
 661
 662   PostData(1, fOutput);
 663 }
 664
 665 //_________________________________________________________________________
 666 void AliAnalysisTaskEmcalJetHadEPpid::ExecOnce()
 667 {
 668   AliAnalysisTaskEmcalJet::ExecOnce();
 669
 670   if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
 671   if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
 672   if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
 673 }
 674
 675 //_________________________________________________________________________
 676 Bool_t AliAnalysisTaskEmcalJetHadEPpid::Run()
 677 { // Main loop called for each event
 678   // TEST TEST TEST TEST for OBJECTS!
 679   if(!fLocalRho){
 680     AliError(Form("Couldn't get fLocalRho object, try to get it from Event based on name\n"));
 681     fLocalRho = GetLocalRhoFromEvent(fLocalRhoName);
 682     if(!fLocalRho) return kTRUE;
 683   }
 684   if(!fTracks){
 685     AliError(Form("No fTracks object!!\n"));
 686     return kTRUE;
 687   }
 688   if(!fJets){
 689     AliError(Form("No fJets object!!\n"));
 690     return kTRUE;
 691   }
 692
 693   // what kind of event do we have: AOD or ESD?
 694   // if we have ESD event, set up ESD object
 695   if(!useAOD){
 696     fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 697     if (!fESD) {
 698       AliError(Form("ERROR: fESD not available\n"));
 699       return kTRUE;
 700     }
 701   }
 702
 703   // if we have AOD event, set up AOD object
 704   if(useAOD){
 705     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 706     if(!fAOD) {
 707       AliError(Form("ERROR: fAOD not available\n"));
 708       return kTRUE;
 709     }
 710   }
 711
 712   // get centrality
 713   Int_t centbin = GetCentBin(fCent);
 714   if (makeQAhistos) fHistCentrality->Fill(fCent); // won't be filled in pp collision (Keep this in mind!)
 715
 716   // for pp analyses we will just use the first centrality bin
 717   if (centbin == -1) centbin = 0;
 718
 719   // apply cut to event on Centrality > 90%
 720   if(fCent>90) return kTRUE;
 721
 722   // get vertex information
 723   Double_t fvertex[3]={0,0,0};
 724   InputEvent()->GetPrimaryVertex()->GetXYZ(fvertex);
 725   Double_t zVtx=fvertex[2];
 726   if (makeQAhistos) fHistZvtx->Fill(zVtx);
 727
 728   // get z-vertex bin
 729   //Int_t zVbin = GetzVertexBin(zVtx);
 730
 731   // apply zVtx cut
 732   if(fabs(zVtx)>10.0) return kTRUE;
 733
 734   // create pointer to list of input event
 735   TList *list = InputEvent()->GetList();
 736   if(!list) {
 737     AliError(Form("ERROR: list not attached\n"));
 738     return kTRUE;
 739   }
 740
 741   // initialize TClonesArray pointers to jets and tracks
 742   TClonesArray *jets = 0;
 743   TClonesArray *tracks = 0;
 744
 745   // get Tracks object
 746   tracks = dynamic_cast<TClonesArray*>(list->FindObject(fTracks));
 747   if (!tracks) {
 748     AliError(Form("Pointer to tracks %s == 0", fTracks->GetName()));
 749     return kTRUE;
 750   } // verify existence of tracks
 751
 752   // get Jets object
 753   jets = dynamic_cast<TClonesArray*>(list->FindObject(fJets));
 754   if(!jets){
 755     AliError(Form("Pointer to jets %s == 0", fJets->GetName()));
 756     return kTRUE;
 757   } // verify existence of jets
 758
 759   // get number of jets and tracks
 760   const Int_t Njets = jets->GetEntries();
 761   const Int_t Ntracks = tracks->GetEntries();
 762   if(Ntracks<1)   return kTRUE;
 763   if(Njets<1)     return kTRUE;
 764
 765   if (makeQAhistos) fHistMult->Fill(Ntracks);  // fill multiplicity distribution
 766
 767   // initialize track parameters
 768   Int_t iTT=-1;
 769   Double_t ptmax=-10;
 770
 771   // loop over tracks - to get hardest track (highest pt)
 772   for (Int_t iTracks = 0; iTracks < Ntracks; iTracks++){
 773     AliVTrack* track = static_cast<AliVTrack*>(tracks->At(iTracks));
 774     if (!track) {
 775       AliError(Form("Couldn't get VTrack track %d\n", iTracks));
 776       continue;
 777     } // verify existence of tracks
 778
 779     // track cuts
 780     if(TMath::Abs(track->Eta())>0.9) continue;
 781     if(track->Pt()<0.15) continue;
 782     //iCount++;
 783     if(track->Pt()>ptmax){
 784       ptmax=track->Pt();             // max pT track
 785       iTT=iTracks;                   // trigger tracks
 786     } // check if Pt>maxpt
 787
 788     if (makeQAhistos) fHistTrackPhi->Fill(track->Phi());
 789     if (makeQAhistos) fHistTrackPt[centbin]->Fill(track->Pt());
 790     if (makeQAhistos) fHistTrackPtallcent->Fill(track->Pt());
 791   } // end of loop over tracks
 792
 793   // get rho from event and fill relative histo's
 794   fRho = GetRhoFromEvent(fRhoName);
 795   fRhoVal = fRho->GetVal();
 796
 797   if (makeQAhistos) {
 798     fHistRhovsdEP[centbin]->Fill(fRhoVal,fEPV0); // Global Rho vs delta Event Plane angle
 799     fHistRhovsCent->Fill(fCent,fRhoVal);        // Global Rho vs Centrality
 800     fHistEP0[centbin]->Fill(fEPV0);
 801     fHistEP0A[centbin]->Fill(fEPV0A);
 802     fHistEP0C[centbin]->Fill(fEPV0C);
 803     fHistEPAvsC[centbin]->Fill(fEPV0A,fEPV0C);
 804   }
 805
 806   // initialize jet parameters
 807   Int_t ijethi=-1;
 808   Double_t highestjetpt=0.0;
 809   Int_t passedTTcut=0;
 810   Int_t NjetAcc = 0;
 811   Double_t leadhadronPT = 0;
 812
 813   // loop over jets in an event - to find highest jet pT and apply some cuts
 814   for (Int_t ijet = 0; ijet < Njets; ijet++){
 815     // get our jets
 816     AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(ijet));
 817     if (!jet) continue;
 818
 819     // apply jet cuts
 820     if ((jet->Phi()<fPhimin)||(jet->Phi()>fPhimax)) continue;
 821     if ((jet->Eta()<fEtamin)||(jet->Eta()>fEtamax)) continue;
 822     if (makeQAhistos) fHistAreavsRawPt[centbin]->Fill(jet->Pt(),jet->Area());
 823     if(!AcceptMyJet(jet)) continue;
 824
 825     NjetAcc++;                     // # of accepted jets
 826
 827     // use this to get total # of jets passing cuts in events!!!!!!!!
 828     if (makeQAhistos) fHistJetPhi->Fill(jet->Phi()); // Jet Phi histogram (filled)
 829
 830     // get highest Pt jet in event
 831     if(highestjetpt<jet->Pt()){
 832       ijethi=ijet;
 833       highestjetpt=jet->Pt();
 834     }
 835   } // end of looping over jets
 836
 837   // accepted jets
 838   fHistNjetvsCent->Fill(fCent,NjetAcc);
 839   Int_t NJETAcc = 0;
 840
 841   // loop over jets in event and make appropriate cuts
 842   for (Int_t ijet = 0; ijet < Njets; ++ijet) {
 843      AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(ijet));
 844      if (!jet) continue;
 845
 846          // (should probably be higher..., but makes a cut on jet pT)
 847      if (jet->Pt()<0.1) continue;
 848      // do we accept jet? apply jet cuts
 849      if (!AcceptMyJet(jet)) continue;
 850
 851      // check on lead jet
 852      Double_t leadjet=0;
 853      if (ijet==ijethi) leadjet=1;
 854
 855      // check on leading hadron pt
 856      leadhadronPT = GetLeadingHadronPt(jet);
 857
 858      // initialize and calculate various parameters: pt, eta, phi, rho, etc...
 859      Double_t jetphi = jet->Phi();      // phi of jet
 860      NJETAcc++;   // # accepted jets
 861      fLocalRhoVal = fLocalRho->GetLocalVal(jetphi, fJetRad); //GetJetRadius(0)); // get local rho value
 862      Double_t jeteta = jet->Eta();     // ETA of jet
 863      Double_t jetPt = -500;
 864      Double_t jetPtGlobal = -500;
 865      Double_t jetPtLocal = -500;            // initialize corr jet pt
 866      jetPt = jet->Pt();
 867      jetPtGlobal = jet->Pt()-jet->Area()*fRhoVal;  // corrected pT of jet from rho value
 868      jetPtLocal = jet->Pt()-jet->Area()*fLocalRhoVal; // corrected pT of jet using Rho modulated for V2 and V3
 869      Double_t dEP = -500;                    // initialize angle between jet and event plane
 870      dEP = RelativeEPJET(jetphi,fEPV0); // angle betweeen jet and event plane
 871
 872      // make histo's
 873      if(makeQAhistos) fHistJetPtvsTrackPt[centbin]->Fill(jetPt,jet->MaxTrackPt());
 874      if(makeQAhistos) fHistRawJetPtvsTrackPt[centbin]->Fill(jetPt,jet->MaxTrackPt());
 875      if(makeQAhistos) fHistJetPtcorrGlRho[centbin]->Fill(jetPtGlobal);
 876      if(makeQAhistos) fHistJetPtvsdEP[centbin]->Fill(jetPt, dEP);
 877      if(makeQAhistos) fHistJetEtaPhiPt[centbin]->Fill(jetPt,jet->Eta(),jet->Phi());
 878      if(makeQAhistos) fHistJetPtArea[centbin]->Fill(jetPt,jet->Area());
 879      if(makeQAhistos) fHistJetEtaPhi->Fill(jet->Eta(),jet->Phi());     // fill jet eta-phi distribution histo
 880      if(makeextraCORRhistos) fHistJetPtNcon[centbin]->Fill(jetPt,jet->GetNumberOfConstituents());
 881      if(makeextraCORRhistos) fHistJetPtNconCh[centbin]->Fill(jetPt,jet->GetNumberOfTracks());
 882      if(makeextraCORRhistos) fHistJetPtNconEm[centbin]->Fill(jetPt,jet->GetNumberOfClusters());
 883      if (makeoldJEThadhistos) fHistJetPt[centbin]->Fill(jet->Pt());  // fill #jets vs pT histo
 884      //fHistDeltaPtvsArea->Fill(jetPt,jet->Area());
 885
 886      // make histo's with BIAS applied
 887      if (jet->MaxTrackPt()>fTrkBias){
 888        if(makeBIAShistos) fHistJetPtvsdEPBias[centbin]->Fill(jetPt, dEP);
 889        if(makeBIAShistos) fHistJetEtaPhiPtBias[centbin]->Fill(jetPt,jet->Eta(),jet->Phi());
 890        if(makeextraCORRhistos) fHistJetPtAreaBias[centbin]->Fill(jetPt,jet->Area());
 891        if(makeextraCORRhistos) fHistJetPtNconBias[centbin]->Fill(jetPt,jet->GetNumberOfConstituents());
 892        if(makeextraCORRhistos) fHistJetPtNconBiasCh[centbin]->Fill(jetPt,jet->GetNumberOfTracks());
 893        if(makeextraCORRhistos) fHistJetPtNconBiasEm[centbin]->Fill(jetPt,jet->GetNumberOfClusters());
 894      }
 895
 896     //if(leadjet && centbin==0){
 897     //  if(makeextraCORRhistos) fHistJetPt[centbin+1]->Fill(jet->Pt());
 898     //}
 899     if ((jet->MaxTrackPt()>fTrkBias) || (jet->MaxClusterPt()>fClusBias)){
 900       if (makeoldJEThadhistos){
 901         fHistJetPtBias[centbin]->Fill(jet->Pt());
 902         //if(leadjet && centbin==0) fHistJetPtBias[centbin+1]->Fill(jet->Pt());
 903       }
 904     }  // end of MaxTrackPt>ftrkBias or maxclusterPt>fclusBias
 905
 906     // do we have trigger tracks
 907     if(iTT>0){
 908       AliVTrack* TT = static_cast<AliVTrack*>(tracks->At(iTT));
 909       if(TMath::Abs(jet->Phi()-TT->Phi()-TMath::Pi())<0.6) passedTTcut=1;
 910       else passedTTcut=0;
 911     } // end of check on iTT > 0
 912
 913     if(passedTTcut) {
 914       if (makeoldJEThadhistos) fHistJetPtTT[centbin]->Fill(jet->Pt());
 915     }
 916
 917     // cut on HIGHEST jet pt in event (15 GeV default)
 918     if (highestjetpt>fJetPtcut) {
 919
 920       // does our max track or cluster pass the bias?
 921       if ((jet->MaxTrackPt()>fTrkBias) || (jet->MaxClusterPt()>fClusBias)){
 922         // set up and fill Jet-Hadron Correction THnSparse
 923         Double_t CorrEntries[3] = {fCent, jet->Pt(), dEP};
 924         fhnCorr->Fill(CorrEntries);    // fill Sparse Histo with Correction entries
 925       }
 926
 927       // loop over all track for an event containing a jet with a pT>fJetPtCut  (15)GeV
 928       for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
 929         AliVTrack* track = static_cast<AliVTrack*>(tracks->At(iTracks));
 930         if (!track) {
 931           AliError(Form("Couldn't get AliVtrack %d\n", iTracks));
 932           continue;
 933         }
 934
 935         // apply track cuts
 936         if(TMath::Abs(track->Eta())>fTrkEta) continue;
 937         if (track->Pt()<0.15) continue;
 938
 939         // calculate and get some track parameters
 940                 Double_t trCharge = -99;
 941         trCharge = track->Charge();
 942         Double_t tracketa=track->Eta();   // eta of track
 943         Double_t deta=tracketa-jeteta;    // dETA between track and jet
 944         //Double_t dR=sqrt(deta*deta+dphijh*dphijh);     // difference of R between jet and hadron track
 945
 946         Int_t ieta = -1;       // initialize deta bin
 947         Int_t iptjet = -1;     // initialize jet pT bin
 948         if (makeoldJEThadhistos)  {
 949                   ieta=GetEtaBin(deta);             // bin of eta
 950               if(ieta<0) continue;              // double check we don't have a negative array index
 951           iptjet=GetpTjetBin(jet->Pt());    // bin of jet pt
 952               if(iptjet<0) continue;                      // double check we don't have a negative array index
 953                 }
 954
 955         // dPHI between jet and hadron
 956         Double_t dphijh = RelativePhi(jet->Phi(), track->Phi()); // angle between jet and hadron
 957
 958         // fill some jet-hadron histo's
 959         if (makeoldJEThadhistos) fHistJetH[centbin][iptjet][ieta]->Fill(dphijh,track->Pt());  // fill jet-hadron dPHI--track pT distribution
 960         if(makeQAhistos) fHistJetHEtaPhi->Fill(deta,dphijh);                          // fill jet-hadron  eta--phi distribution
 961             fHistJetHaddPHI->Fill(dphijh);
 962         if(passedTTcut){
 963           if (makeoldJEThadhistos) fHistJetHTT[centbin][iptjet][ieta]->Fill(dphijh,track->Pt());
 964         }
 965
 966         // does our max track or cluster pass the bias?
 967         if ((jet->MaxTrackPt()>fTrkBias) || (jet->MaxClusterPt()>fClusBias)){
 968           // set up and fill Jet-Hadron THnSparse
 969           Double_t triggerEntries[9] = {fCent, jet->Pt(), track->Pt(), deta, dphijh, dEP, zVtx, trCharge, leadjet};
 970           fhnJH->Fill(triggerEntries);    // fill Sparse Histo with trigger entries
 971
 972           // fill histo's
 973           if(makeQAhistos) fHistSEphieta->Fill(dphijh, deta); // single event distribution
 974           if (makeoldJEThadhistos) fHistJetHBias[centbin][iptjet][ieta]->Fill(dphijh,track->Pt());
 975
 976           if (makeBIAShistos) {
 977             fHistJetHaddPhiBias->Fill(dphijh);
 978
 979             // in plane and out of plane histo's
 980             if( dEP>0 && dEP<=(TMath::Pi()/6) ){
 981               // we are IN plane
 982               fHistJetHaddPhiINBias->Fill(dphijh);
 983             }else if( dEP>(TMath::Pi()/3) && dEP<=(TMath::Pi()/2) ){
 984               // we are OUT of PLANE
 985               fHistJetHaddPhiOUTBias->Fill(dphijh);
 986             }else if( dEP>(TMath::Pi()/6) && dEP<=(TMath::Pi()/3) ){
 987               // we are in middle of plane
 988               fHistJetHaddPhiMIDBias->Fill(dphijh);
 989             }
 990                   } // BIAS histos switch
 991         } // end of check maxtrackpt>ftrackbias or maxclusterpt>fclustbias
 992
 993         // **************************************************************************************************************
 994         // *********************************** PID **********************************************************************
 995         // **************************************************************************************************************
 996         if(doPIDtrackBIAS){
 997           //if(ptmax < fTrkBias) continue;    // force PID to happen when max track pt > 5.0 GeV
 998           if(leadhadronPT < fTrkBias) continue; // force PID to happen when lead hadron pt > 5.0 GeV
 999         }
1000
1001         // some variables for PID
1002             Double_t eta = -999;
1003         Double_t pt = -999;
1004         Double_t dEdx = -999;
1005         Double_t ITSsig = -999;
1006         Double_t TOFsig = -999;
1007         Double_t charge = -999;
1008
1009         // nSigma of particles in TPC, TOF, and ITS
1010         Double_t nSigmaPion_TPC, nSigmaProton_TPC, nSigmaKaon_TPC;
1011         Double_t nSigmaPion_TOF, nSigmaProton_TOF, nSigmaKaon_TOF;
1012         Double_t nSigmaPion_ITS, nSigmaProton_ITS, nSigmaKaon_ITS;
1013
1014         if(doPID){
1015           // get parameters of track
1016           charge = track->Charge();    // charge of track
1017           eta    = track->Eta();       // ETA of track
1018           pt     = track->Pt();        // pT of track
1019
1020           // extra attempt
1021           AliVEvent *vevent=InputEvent();
1022           if (!vevent||!fPIDResponse) return kTRUE; // just return, maybe put at beginning
1023
1024           // get PID parameters, first check if AOD/ESD
1025               if (!useAOD) {
1026             AliESDtrack *trackESD = fESD->GetTrack(iTracks);
1027
1028             // get detector signals
1029             dEdx = trackESD->GetTPCsignal();
1030             ITSsig = trackESD->GetITSsignal();
1031             TOFsig = trackESD->GetTOFsignal();
1032
1033             // TPC nSigma's
1034             nSigmaPion_TPC = fPIDResponse->NumberOfSigmasTPC(trackESD,AliPID::kPion);
1035             nSigmaKaon_TPC = fPIDResponse->NumberOfSigmasTPC(trackESD,AliPID::kKaon);
1036             nSigmaProton_TPC = fPIDResponse->NumberOfSigmasTPC(trackESD,AliPID::kProton);
1037
1038             // TOF nSigma's
1039             nSigmaPion_TOF = fPIDResponse->NumberOfSigmasTOF(trackESD,AliPID::kPion);
1040             nSigmaKaon_TOF = fPIDResponse->NumberOfSigmasTOF(trackESD,AliPID::kKaon);
1041             nSigmaProton_TOF = fPIDResponse->NumberOfSigmasTOF(trackESD,AliPID::kProton);
1042
1043             // ITS nSigma's
1044             nSigmaPion_ITS = fPIDResponse->NumberOfSigmasITS(trackESD,AliPID::kPion);
1045             nSigmaKaon_ITS = fPIDResponse->NumberOfSigmasITS(trackESD,AliPID::kKaon);
1046             nSigmaProton_ITS = fPIDResponse->NumberOfSigmasITS(trackESD,AliPID::kProton);
1047               } // end of ESD pid
1048
1049           if (useAOD) {
1050                 AliAODTrack *trackAOD = fAOD->GetTrack(iTracks);
1051
1052             // get detector signals
1053             dEdx = trackAOD->GetTPCsignal();
1054             ITSsig = trackAOD->GetITSsignal();
1055             TOFsig = trackAOD->GetTOFsignal();
1056
1057             // TPC nSigma's
1058             nSigmaPion_TPC = fPIDResponse->NumberOfSigmasTPC(trackAOD,AliPID::kPion);
1059             nSigmaKaon_TPC = fPIDResponse->NumberOfSigmasTPC(trackAOD,AliPID::kKaon);
1060             nSigmaProton_TPC = fPIDResponse->NumberOfSigmasTPC(trackAOD,AliPID::kProton);
1061
1062             // TOF nSigma's
1063             nSigmaPion_TOF = fPIDResponse->NumberOfSigmasTOF(trackAOD,AliPID::kPion);
1064             nSigmaKaon_TOF = fPIDResponse->NumberOfSigmasTOF(trackAOD,AliPID::kKaon);
1065             nSigmaProton_TOF = fPIDResponse->NumberOfSigmasTOF(trackAOD,AliPID::kProton);
1066
1067             // ITS nSigma's
1068             nSigmaPion_ITS = fPIDResponse->NumberOfSigmasITS(trackAOD,AliPID::kPion);
1069             nSigmaKaon_ITS = fPIDResponse->NumberOfSigmasITS(trackAOD,AliPID::kKaon);
1070             nSigmaProton_ITS = fPIDResponse->NumberOfSigmasITS(trackAOD,AliPID::kProton);
1071               } // end of AOD pid
1072
1073           // fill detector signal histograms
1074           if (makeQAhistos) fHistTPCdEdX->Fill(pt, dEdx);
1075           if (makeQAhistos) fHistITSsignal->Fill(pt, ITSsig);
1076           //if (makeQAhistos) fHistTOFsignal->Fill(pt, TOFsig);
1077
1078           // Tests to PID pions, kaons, and protons,          (default is undentified tracks)
1079           Double_t nPIDtpc = 0;
1080           Double_t nPIDits = 0;
1081           Double_t nPIDtof = 0;
1082           Double_t nPID = -99;
1083
1084           // check track has pT < 0.900 GeV  - use TPC pid
1085           if (pt<0.900 && dEdx>0) {
1086                 nPIDtpc = 4;
1087             nPID = 0.5;
1088
1089             // PION check - TPC
1090             if (TMath::Abs(nSigmaPion_TPC)<2 && TMath::Abs(nSigmaKaon_TPC)>2 && TMath::Abs(nSigmaProton_TPC)>2 ){
1091               isPItpc = kTRUE;
1092               nPIDtpc = 1;
1093               nPID=1.5;
1094             }else isPItpc = kFALSE;
1095
1096             // KAON check - TPC
1097             if (TMath::Abs(nSigmaKaon_TPC)<2 && TMath::Abs(nSigmaPion_TPC)>3 && TMath::Abs(nSigmaProton_TPC)>2 ){
1098               isKtpc = kTRUE;
1099               nPIDtpc = 2;
1100               nPID=2.5;
1101             }else isKtpc = kFALSE;
1102
1103             // PROTON check - TPC
1104             if (TMath::Abs(nSigmaProton_TPC)<2 && TMath::Abs(nSigmaPion_TPC)>3 && TMath::Abs(nSigmaKaon_TPC)>2 ){
1105               isPtpc = kTRUE;
1106               nPIDtpc = 3;
1107               nPID=3.5;
1108             }else isPtpc = kFALSE;
1109           }  // cut on track pT for TPC
1110
1111           // check track has pT < 0.500 GeV - use ITS pid
1112           if (pt<0.500 && ITSsig>0) {
1113             nPIDits = 4;
1114             nPID = 4.5;
1115
1116             // PION check - ITS
1117             if (TMath::Abs(nSigmaPion_ITS)<2 && TMath::Abs(nSigmaKaon_ITS)>2 && TMath::Abs(nSigmaProton_ITS)>2 ){
1118               isPIits = kTRUE;
1119               nPIDits = 1;
1120                   nPID=5.5;
1121             }else isPIits = kFALSE;
1122
1123             // KAON check - ITS
1124             if (TMath::Abs(nSigmaKaon_ITS)<2 && TMath::Abs(nSigmaPion_ITS)>3 && TMath::Abs(nSigmaProton_ITS)>2 ){
1125               isKits = kTRUE;
1126               nPIDits = 2;
1127               nPID=6.5;
1128             }else isKits = kFALSE;
1129
1130             // PROTON check - ITS
1131             if (TMath::Abs(nSigmaProton_ITS)<2 && TMath::Abs(nSigmaPion_ITS)>3 && TMath::Abs(nSigmaKaon_ITS)>2 ){
1132               isPits = kTRUE;
1133               nPIDits = 3;
1134                   nPID=7.5;
1135             }else isPits = kFALSE;
1136           }  // cut on track pT for ITS
1137
1138           // check track has 0.900 GeV < pT < 2.500 GeV - use TOF pid
1139           if (pt>0.900 && pt<2.500 && TOFsig>0) {
1140                 nPIDtof = 4;
1141             nPID = 8.5;
1142
1143             // PION check - TOF
1144             if (TMath::Abs(nSigmaPion_TOF)<2 && TMath::Abs(nSigmaKaon_TOF)>2 && TMath::Abs(nSigmaProton_TOF)>2 ){
1145               isPItof = kTRUE;
1146               nPIDtof = 1;
1147               nPID=9.5;
1148             }else isPItof = kFALSE;
1149
1150             // KAON check - TOF
1151             if (TMath::Abs(nSigmaKaon_TOF)<2 && TMath::Abs(nSigmaPion_TOF)>3 && TMath::Abs(nSigmaProton_TOF)>2 ){
1152               isKtof = kTRUE;
1153               nPIDtof = 2;
1154               nPID=10.5;
1155             }else isKtof = kFALSE;
1156
1157             // PROTON check - TOF
1158             if (TMath::Abs(nSigmaProton_TOF)<2 && TMath::Abs(nSigmaPion_TOF)>3 && TMath::Abs(nSigmaKaon_TOF)>2 ){
1159               isPtof = kTRUE;
1160               nPIDtof = 3;
1161               nPID=11.5;
1162             }else isPtof = kFALSE;
1163           }  // cut on track pT for TOF
1164
1165           if (nPID == -99) nPID = 13.5;
1166               fHistPID->Fill(nPID);
1167
1168           // PID sparse getting filled
1169           Double_t pid_Entries[23] = {fCent,jetPtLocal,pt,charge,eta,deta,dphijh,leadjet,zVtx,dEP,jetPt,
1170                                       nSigmaPion_TPC, nSigmaProton_TPC, nSigmaKaon_TPC,  //nSig in TPC
1171                                       nSigmaPion_ITS, nSigmaProton_ITS, nSigmaKaon_ITS,  //nSig in ITS
1172                                       nSigmaPion_TOF, nSigmaProton_TOF, nSigmaKaon_TOF,  //nSig in TOF
1173                                                                   nPIDtpc, nPIDits, nPIDtof                          //PID label for each detector
1174                                       }; //array for PID sparse
1175           fhnPID->Fill(pid_Entries);   // fill Sparse histo of PID tracks
1176         } // end of doPID check
1177
1178             // get track pt bin
1179         Int_t itrackpt = -500;              // initialize track pT bin
1180         itrackpt = GetpTtrackBin(track->Pt());
1181
1182             // all tracks: jet hadron correlations in hadron pt bins
1183         if(makeextraCORRhistos) fHistJetHadbindPhi[itrackpt]->Fill(dphijh);
1184
1185         // in plane and out of plane jet-hadron histo's
1186         if( dEP>0 && dEP<=(TMath::Pi()/6) ){
1187           // we are IN plane
1188           if(makeextraCORRhistos) fHistJetHaddPhiINcent[centbin]->Fill(dphijh);
1189           fHistJetHaddPhiIN->Fill(dphijh);
1190           if(makeextraCORRhistos) fHistJetHadbindPhiIN[itrackpt]->Fill(dphijh);
1191           //fHistJetHaddPhiPtcentbinIN[itrackpt][centbin]->Fill(dphijh);
1192         }else if( dEP>(TMath::Pi()/3) && dEP<=(TMath::Pi()/2) ){
1193           // we are OUT of PLANE
1194           if(makeextraCORRhistos) fHistJetHaddPhiOUTcent[centbin]->Fill(dphijh);
1195           fHistJetHaddPhiOUT->Fill(dphijh);
1196           if(makeextraCORRhistos) fHistJetHadbindPhiOUT[itrackpt]->Fill(dphijh);
1197           //fHistJetHaddPhiPtcentbinOUT[itrackpt][centbin]->Fill(dphijh);
1198         }else if( dEP>(TMath::Pi()/6) && dEP<=(TMath::Pi()/3) ){
1199           // we are in the middle of plane
1200           if(makeextraCORRhistos) fHistJetHaddPhiMIDcent[centbin]->Fill(dphijh);
1201           fHistJetHaddPhiMID->Fill(dphijh);
1202           if(makeextraCORRhistos) fHistJetHadbindPhiMID[itrackpt]->Fill(dphijh);
1203         }
1204       } // loop over tracks found in event with highest JET pT > 10.0 GeV (change)
1205     } // jet pt cut
1206   } // jet loop
1207
1208 // ***************************************************************************************************************
1209 // ******************************** Event MIXING *****************************************************************
1210   TObjArray* tracksClone = CloneAndReduceTrackList(tracks); // TEST
1211
1212   //Prepare to do event mixing
1213   if(fDoEventMixing>0){
1214     // event mixing
1215
1216     // 1. First get an event pool corresponding in mult (cent) and
1217     //    zvertex to the current event. Once initialized, the pool
1218     //    should contain nMix (reduced) events. This routine does not
1219     //    pre-scan the chain. The first several events of every chain
1220     //    will be skipped until the needed pools are filled to the
1221     //    specified depth. If the pool categories are not too rare, this
1222     //    should not be a problem. If they are rare, you could lose
1223     //    statistics.
1224
1225     // 2. Collect the whole pool's content of tracks into one TObjArray
1226     //    (bgTracks), which is effectively a single background super-event.
1227
1228     // 3. The reduced and bgTracks arrays must both be passed into
1229     //    FillCorrelations(). Also nMix should be passed in, so a weight
1230     //    of 1./nMix can be applied.
1231
1232     // mix jets from triggered events with tracks from MB events
1233     // get the trigger bit
1234     // need to change trigger bits between different runs
1235 //T    UInt_t trigger = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
1236 //T    if (trigger==0)  return kTRUE; // return
1237
1238     //Double_t Ntrks=(Double_t)Ntracks*1.0;
1239     //cout<<"Test.. Ntrks: "<<fPoolMgr->GetEventPool(Ntrks);
1240
1241     AliEventPool* pool = fPoolMgr->GetEventPool(fCent, zVtx); // for PbPb? fcent
1242     //AliEventPool* pool = fPoolMgr->GetEventPool(Ntrks, zVtx); // for pp
1243
1244     if (!pool){
1245       AliFatal(Form("No pool found for centrality = %f, zVtx = %f", fCent, zVtx));
1246       //AliFatal(Form("No pool found for multiplicity = %f, zVtx = %f", Ntrks, zVtx));
1247       return kTRUE;
1248     }
1249
1250     // use only jets from EMCal-triggered events (for lhc11a use AliVEvent::kEMC1)
1251 ///    if (trigger & AliVEvent::kEMC1) {
1252 //T    if (trigger & AliVEvent::kEMCEJE) {  // TEST
1253       //check for a trigger jet
1254       // fmixingtrack/10 ??
1255       if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 || pool->GetCurrentNEvents() >= 5) {
1256         // loop over jets (passing cuts?)
1257         for (Int_t ijet = 0; ijet < Njets; ijet++) {
1258           Double_t leadjet=0;
1259           if (ijet==ijethi) leadjet=1;
1260
1261           // get jet object
1262           AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(ijet));
1263                   if (!jet) continue;
1264
1265                   // (should probably be higher..., but makes a cut on jet pT)
1266           if (jet->Pt()<0.1) continue;
1267           if (!AcceptMyJet(jet)) continue;
1268
1269           Int_t nMix = pool->GetCurrentNEvents();  // how many particles in pool to mix
1270
1271           // Fill for biased jet triggers only
1272           if ((jet->MaxTrackPt()>fTrkBias) || (jet->MaxClusterPt()>fClusBias)) {
1273             // Fill mixed-event histos here
1274             for (Int_t jMix=0; jMix<nMix; jMix++) {
1275                 TObjArray* bgTracks = pool->GetEvent(jMix);
1276                 const Int_t Nbgtrks = bgTracks->GetEntries();
1277                 for(Int_t ibg=0; ibg<Nbgtrks; ibg++) {
1278                   AliPicoTrack *part = static_cast<AliPicoTrack*>(bgTracks->At(ibg));
1279                   if(!part) continue;
1280                   if(TMath::Abs(part->Eta())>0.9) continue;
1281                   if(part->Pt()<0.15) continue;
1282
1283                   Double_t DEta = part->Eta()-jet->Eta();                // difference in eta
1284                   Double_t DPhi = RelativePhi(jet->Phi(),part->Phi());   // difference in phi
1285                   Double_t dEP = RelativeEPJET(jet->Phi(),fEPV0);            // difference between jet and EP
1286                           Double_t mixcharge = part->Charge();
1287                   //Double_t DR=TMath::Sqrt(DPhi*DPhi+DEta*DEta);      // difference in R
1288
1289                   Double_t triggerEntries[10] = {fCent,jet->Pt(),part->Pt(),DEta,DPhi,dEP,zVtx, mixcharge, leadjet}; //array for ME sparse
1290                   fhnMixedEvents->Fill(triggerEntries,1./nMix);   // fill Sparse histo of mixed events
1291
1292                   if(makeextraCORRhistos) fHistMEphieta->Fill(DPhi,DEta, 1./nMix);
1293                 } // end of background track loop
1294              } // end of filling mixed-event histo's
1295           } // end of check for biased jet triggers
1296         } // end of jet loop
1297       } // end of check for triggered jet
1298 //    } //end EMC triggered loop
1299
1300     // use only tracks from MB events (for lhc11a use AliVEvent::kMB)
1301 ///    if (trigger & AliVEvent::kMB) {
1302 //T    if (trigger & AliVEvent::kAnyINT){ // test
1303       // create a list of reduced objects. This speeds up processing and reduces memory consumption for the event pool
1304 //T      TObjArray* tracksClone = CloneAndReduceTrackList(tracks);
1305
1306       // update pool if jet in event or not
1307       pool->UpdatePool(tracksClone);
1308 ///    } // check on track from MB events
1309   } // end of event mixing
1310
1311   // print some stats on the event
1312   event++;
1313
1314   if (doComments) {
1315     cout<<"Event #: "<<event<<"     Jet Radius: "<<fJetRad<<"     Constituent Pt Cut: "<<fConstituentCut<<endl;
1316     cout<<"# of jets: "<<Njets<<"      Highest jet pt: "<<highestjetpt<<"     leading hadron pt: "<<leadhadronPT<<endl;
1317     cout<<"# tracks: "<<Ntracks<<"      Highest track pt: "<<ptmax<<endl;
1318     cout<<" =============================================== "<<endl;
1319   }
1320
1321   return kTRUE;  // used when the function is of type bool
1322 }  // end of RUN
1323
1324 //________________________________________________________________________
1325 Int_t AliAnalysisTaskEmcalJetHadEPpid::GetCentBin(Double_t cent) const
1326 {  // Get centrality bin.
1327   Int_t centbin = -1;
1328   if (cent>=0 && cent<10)       centbin = 0;
1329   else if (cent>=10 && cent<20) centbin = 1;
1330   else if (cent>=20 && cent<30) centbin = 2;
1331   else if (cent>=30 && cent<40) centbin = 3;
1332   else if (cent>=40 && cent<50) centbin = 4;
1333   else if (cent>=50 && cent<90) centbin = 5;
1334
1335   return centbin;
1336 }
1337
1338 //________________________________________________________________________
1339 Double_t AliAnalysisTaskEmcalJetHadEPpid::RelativePhi(Double_t mphi,Double_t vphi) const
1340 { // function to calculate relative PHI
1341   double dphi = mphi-vphi;
1342
1343   // set dphi to operate on adjusted scale
1344   if(dphi<-0.5*TMath::Pi()) dphi+=2.*TMath::Pi();
1345   if(dphi>3./2.*TMath::Pi()) dphi-=2.*TMath::Pi();
1346
1347   // test
1348   if( dphi < -1.*TMath::Pi()/2 || dphi > 3.*TMath::Pi()/2 )
1349     AliWarning(Form("%s: dPHI not in range [-0.5*Pi, 1.5*Pi]!", GetName()));
1350
1351   return dphi; // dphi in [-0.5Pi, 1.5Pi]
1352 }
1353
1354
1355
1356 //_________________________________________________________________________
1357 Double_t AliAnalysisTaskEmcalJetHadEPpid:: RelativeEPJET(Double_t jetAng, Double_t EPAng) const
1358 { // function to calculate angle between jet and EP in the 1st quadrant (0,Pi/2)
1359   Double_t dphi = (EPAng - jetAng);
1360
1361   // ran into trouble with a few dEP<-Pi so trying this...
1362   if( dphi<-1*TMath::Pi() ){
1363     dphi = dphi + 1*TMath::Pi();
1364   } // this assumes we are doing full jets currently
1365
1366   if( (dphi>0) && (dphi<1*TMath::Pi()/2) ){
1367     // Do nothing! we are in quadrant 1
1368   }else if( (dphi>1*TMath::Pi()/2) && (dphi<1*TMath::Pi()) ){
1369     dphi = 1*TMath::Pi() - dphi;
1370   }else if( (dphi<0) && (dphi>-1*TMath::Pi()/2) ){
1371     dphi = fabs(dphi);
1372   }else if( (dphi<-1*TMath::Pi()/2) && (dphi>-1*TMath::Pi()) ){
1373     dphi = dphi + 1*TMath::Pi();
1374   }
1375
1376   // test
1377   if( dphi < 0 || dphi > TMath::Pi()/2 )
1378     AliWarning(Form("%s: dPHI not in range [0, 0.5*Pi]!", GetName()));
1379
1380   return dphi;   // dphi in [0, Pi/2]
1381 }
1382
1383 //Int_t ieta=GetEtaBin(deta);
1384 //________________________________________________________________________
1385 Int_t AliAnalysisTaskEmcalJetHadEPpid::GetEtaBin(Double_t eta) const
1386 {
1387   // Get eta bin for histos.
1388   Int_t etabin = -1;
1389   if (TMath::Abs(eta)<=0.4)                                             etabin = 0;
1390   else if (TMath::Abs(eta)>0.4 && TMath::Abs(eta)<0.8)  etabin = 1;
1391   else if (TMath::Abs(eta)>=0.8)                                    etabin = 2;
1392
1393   return etabin;
1394 } // end of get-eta-bin
1395
1396 //________________________________________________________________________
1397 Int_t AliAnalysisTaskEmcalJetHadEPpid::GetpTjetBin(Double_t pt) const
1398 {
1399   // Get jet pt  bin for histos.
1400   Int_t ptbin = -1;
1401   if (pt>=15 && pt<20)          ptbin = 0;
1402   else if (pt>=20 && pt<25)     ptbin = 1;
1403   else if (pt>=25 && pt<40)     ptbin = 2;
1404   else if (pt>=40 && pt<60)     ptbin = 3;
1405   else if (pt>=60)              ptbin = 4;
1406
1407   return ptbin;
1408 } // end of get-jet-pt-bin
1409
1410 //________________________________________________________________________
1411 Int_t AliAnalysisTaskEmcalJetHadEPpid::GetpTtrackBin(Double_t pt) const
1412 {
1413   // May need to update bins for future runs... (testing locally)
1414
1415   // Get track pt bin for histos.
1416   Int_t ptbin = -1;
1417   if (pt < 0.5)                 ptbin = 0;
1418   else if (pt>=0.5 && pt<1.0)   ptbin = 1;
1419   else if (pt>=1.0 && pt<1.5)   ptbin = 2;
1420   else if (pt>=1.5 && pt<2.0)   ptbin = 3;
1421   else if (pt>=2.0 && pt<2.5)   ptbin = 4;
1422   else if (pt>=2.5 && pt<3.0)   ptbin = 5;
1423   else if (pt>=3.0 && pt<4.0)   ptbin = 6;
1424   else if (pt>=4.0 && pt<5.0)   ptbin = 7;
1425   else if (pt>=5.0)             ptbin = 8;
1426
1427   return ptbin;
1428 } // end of get-jet-pt-bin
1429
1430
1431 //___________________________________________________________________________
1432 Int_t AliAnalysisTaskEmcalJetHadEPpid::GetzVertexBin(Double_t zVtx) const
1433 {
1434   // get z-vertex bin for histo.
1435   int zVbin= -1;
1436   if (zVtx>=-10 && zVtx<-8)         zVbin = 0;
1437   else if (zVtx>=-8 && zVtx<-6) zVbin = 1;
1438   else if (zVtx>=-6 && zVtx<-4) zVbin = 2;
1439   else if (zVtx>=-4 && zVtx<-2) zVbin = 3;
1440   else if (zVtx>=-2 && zVtx<0)  zVbin = 4;
1441   else if (zVtx>=0 && zVtx<2)   zVbin = 5;
1442   else if (zVtx>=2 && zVtx<4)   zVbin = 6;
1443   else if (zVtx>=4 && zVtx<6)   zVbin = 7;
1444   else if (zVtx>=6 && zVtx<8)   zVbin = 8;
1445   else if (zVtx>=8 && zVtx<10)  zVbin = 9;
1446   else zVbin = 10;
1447
1448   return zVbin;
1449 } // end of get z-vertex bin
1450
1451 //______________________________________________________________________
1452 THnSparse* AliAnalysisTaskEmcalJetHadEPpid::NewTHnSparseD(const char* name, UInt_t entries)
1453 {
1454    // generate new THnSparseD, axes are defined in GetDimParams()
1455    Int_t count = 0;
1456    UInt_t tmp = entries;
1457    while(tmp!=0){
1458       count++;
1459       tmp = tmp &~ -tmp;  // clear lowest bit
1460    }
1461
1462    TString hnTitle(name);
1463    const Int_t dim = count;
1464    Int_t nbins[dim];
1465    Double_t xmin[dim];
1466    Double_t xmax[dim];
1467
1468    Int_t i=0;
1469    Int_t c=0;
1470    while(c<dim && i<32){
1471       if(entries&(1<<i)){
1472
1473          TString label("");
1474          GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
1475          hnTitle += Form(";%s",label.Data());
1476          c++;
1477       }
1478
1479       i++;
1480    }
1481    hnTitle += ";";
1482
1483    return new THnSparseD(name, hnTitle.Data(), dim, nbins, xmin, xmax);
1484 } // end of NewTHnSparseD
1485
1486 void AliAnalysisTaskEmcalJetHadEPpid::GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
1487 {
1488    // stores label and binning of axis for THnSparse
1489    const Double_t pi = TMath::Pi();
1490
1491    switch(iEntry){
1492
1493    case 0:
1494       label = "V0 centrality (%)";
1495       nbins = 10;
1496       xmin = 0.;
1497       xmax = 100.;
1498       break;
1499
1500    case 1:
1501       label = "Jet p_{T}";
1502       nbins = 216;
1503       xmin = 0.;
1504       xmax = 216.;
1505       break;
1506
1507    case 2:
1508       label = "Track p_{T}";
1509       nbins = 750;
1510       xmin = 0.;
1511       xmax = 75.;
1512       break;
1513
1514     case 3:
1515       label = "Relative Eta";
1516       nbins = 48;
1517       xmin = -1.6;
1518       xmax = 1.6;
1519       break;
1520
1521    case 4:
1522       label = "Relative Phi";
1523       nbins = 72;
1524       xmin = -0.5*pi;
1525       xmax = 1.5*pi;
1526       break;
1527
1528   case 5:
1529       label = "Relative angle of Jet and Reaction Plane";
1530       nbins = 72;
1531       xmin = -0.5*pi;
1532       xmax = 1.5*pi;
1533       break;
1534
1535   case 6:
1536       label = "z-vertex";
1537       nbins = 10;
1538       xmin = -10;
1539       xmax =  10;
1540       break;
1541
1542   case 7:
1543       label = "track charge";
1544       nbins = 3;
1545       xmin = -1.5;
1546       xmax = 1.5;
1547       break;
1548
1549   case 8:
1550       label = "leading jet";
1551       nbins = 3;
1552       xmin = -0.5;
1553       xmax = 2.5;
1554       break;
1555
1556   case 9: // need to update
1557       label = "leading track";
1558       nbins = 10;
1559       xmin = 0;
1560       xmax = 50;
1561       break;
1562
1563    } // end of switch
1564 } // end of getting dim-params
1565
1566 //_________________________________________________
1567 // From CF event mixing code PhiCorrelations
1568 TObjArray* AliAnalysisTaskEmcalJetHadEPpid::CloneAndReduceTrackList(TObjArray* tracks)
1569 {
1570   // clones a track list by using AliPicoTrack which uses much less memory (used for event mixing)
1571   TObjArray* tracksClone = new TObjArray;
1572   tracksClone->SetOwner(kTRUE);
1573
1574   /*
1575   Int_t nTrax = fESD->GetNumberOfTracks();
1576   //cout << "nTrax " << nTrax <<endl;
1577
1578   for (Int_t i = 0; i < nTrax; ++i)
1579     {
1580       AliESDtrack* esdtrack = fESD->GetTrack(i);
1581       if (!esdtrack)
1582         {
1583           AliError(Form("Couldn't get ESD track %d\n", i));
1584           continue;
1585         }
1586
1587       AliESDtrack *particle = GetAcceptTrack(esdtrack);
1588       if(!particle) continue;
1589
1590       cout << "RM Hybrid track : " << i << "  " << particle->Pt() << endl;
1591   */
1592
1593   //cout << "nEntries " << tracks->GetEntriesFast() <<endl;
1594   for (Int_t i=0; i<tracks->GetEntriesFast(); i++) {
1595     AliVParticle* particle = (AliVParticle*) tracks->At(i);
1596     if(TMath::Abs(particle->Eta())>fTrkEta) continue;
1597     if(particle->Pt()<0.15)continue;
1598
1599 /*
1600 // DON'T USE
1601     Double_t trackpt=particle->Pt();   // track pT
1602
1603     Int_t trklabel=-1;
1604     trklabel=particle->GetLabel();
1605     //cout << "TRACK_LABEL: " << particle->GetLabel()<<endl;
1606
1607     Int_t hadbin=-1;
1608     if(trackpt<0.5) hadbin=0;
1609     else if(trackpt<1) hadbin=1;
1610     else if(trackpt<2) hadbin=2;
1611     else if(trackpt<3) hadbin=3;
1612     else if(trackpt<5) hadbin=4;
1613     else if(trackpt<8) hadbin=5;
1614     else if(trackpt<20) hadbin=6;
1615 // end of DON'T USE
1616
1617 //feb10 comment out
1618     if(hadbin>-1 && trklabel>-1 && trklabel <3) fHistTrackEtaPhi[trklabel][hadbin]->Fill(particle->Eta(),particle->Phi());
1619     if(hadbin>-1) fHistTrackEtaPhi[3][hadbin]->Fill(particle->Eta(),particle->Phi());
1620
1621     if(hadbin>-1) fHistTrackEtaPhi[hadbin]->Fill(particle->Eta(),particle->Phi());  // TEST
1622 */
1623
1624     tracksClone->Add(new AliPicoTrack(particle->Pt(), particle->Eta(), particle->Phi(), particle->Charge(), 0, 0, 0, 0));
1625   } // end of looping through tracks
1626
1627   return tracksClone;
1628 }
1629
1630 //____________________________________________________________________________________________
1631 THnSparse* AliAnalysisTaskEmcalJetHadEPpid::NewTHnSparseDPID(const char* name, UInt_t entries)
1632 {
1633    // generate new THnSparseD PID, axes are defined in GetDimParams()
1634    Int_t count = 0;
1635    UInt_t tmp = entries;
1636    while(tmp!=0){
1637       count++;
1638       tmp = tmp &~ -tmp;  // clear lowest bit
1639    }
1640
1641    TString hnTitle(name);
1642    const Int_t dim = count;
1643    Int_t nbins[dim];
1644    Double_t xmin[dim];
1645    Double_t xmax[dim];
1646
1647    Int_t i=0;
1648    Int_t c=0;
1649    while(c<dim && i<32){
1650       if(entries&(1<<i)){
1651
1652          TString label("");
1653          GetDimParamsPID(i, label, nbins[c], xmin[c], xmax[c]);
1654          hnTitle += Form(";%s",label.Data());
1655          c++;
1656       }
1657
1658       i++;
1659    }
1660    hnTitle += ";";
1661
1662    return new THnSparseD(name, hnTitle.Data(), dim, nbins, xmin, xmax);
1663 } // end of NewTHnSparseD PID
1664
1665 //________________________________________________________________________________
1666 void AliAnalysisTaskEmcalJetHadEPpid::GetDimParamsPID(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
1667 {
1668    // stores label and binning of axis for THnSparse
1669    const Double_t pi = TMath::Pi();
1670
1671    switch(iEntry){
1672
1673    case 0:
1674       label = "V0 centrality (%)";
1675       nbins = 10;
1676       xmin = 0.;
1677       xmax = 100.;
1678       break;
1679
1680    case 1:
1681       label = "Corrected Jet p_{T}";
1682       nbins = 216;
1683       xmin = -50.;
1684       xmax = 166.;
1685       break;
1686
1687    case 2:
1688       label = "Track p_{T}";
1689       nbins = 750;
1690       xmin = 0.;
1691       xmax = 75.;
1692       break;
1693
1694    case 3:
1695       label = "Charge of Track";
1696       nbins = 3;
1697       xmin = -1.5;
1698       xmax = 1.5;
1699       break;
1700
1701    case 4:
1702       label = "Track Eta";
1703       nbins = 48;
1704       xmin = -1.2;
1705       xmax = 1.2;
1706       break;
1707
1708    case 5:
1709       label = "Relative Eta of Track and Jet";
1710       nbins = 48;
1711       xmin = -1.6;
1712       xmax = 1.6;
1713       break;
1714
1715    case 6:
1716       label = "Relative Phi of Track and Jet";
1717       nbins = 72;
1718       xmin = -0.5*pi;
1719       xmax = 1.5*pi;
1720       break;
1721
1722    case 7:
1723       label = "leading jet";
1724       nbins = 3;
1725       xmin = -.5;
1726       xmax = 2.5;
1727       break;
1728
1729    case 8:
1730       label = "Z-vertex";
1731       nbins = 10;
1732       xmin = -10.;
1733       xmax = 10.;
1734       break;
1735
1736    case 9:
1737       label = "Relative angle: Jet and Reaction Plane";
1738       nbins = 48;
1739       xmin = 0.;
1740       xmax = 0.5*pi;
1741       break;
1742
1743    case 10:
1744       label = "Jet Pt";
1745       nbins = 216;
1746       xmin = 0.;
1747       xmax = 216.;
1748       break;
1749
1750    case 11:
1751       label = "N-Sigma of pions in TPC";
1752       nbins = 200;
1753       xmin = -10.0;
1754       xmax = 10.0;
1755       break;
1756
1757    case 12:
1758       label = "N-Sigma of protons in TPC";
1759       nbins = 200;
1760       xmin = -10.;
1761       xmax = 10.;
1762       break;
1763
1764    case 13:
1765       label = "N-Sigma of kaons in TPC";
1766       nbins = 200;
1767       xmin = -10.;
1768       xmax = 10.;
1769       break;
1770
1771    case 14:
1772       label = "N-Sigma of pions in ITS";
1773       nbins = 200;
1774       xmin = -10.0;
1775       xmax = 10.0;
1776       break;
1777
1778    case 15:
1779       label = "N-Sigma of protons in ITS";
1780       nbins = 200;
1781       xmin = -10.;
1782       xmax = 10.;
1783       break;
1784
1785    case 16:
1786       label = "N-Sigma of kaons in ITS";
1787       nbins = 200;
1788       xmin = -10.;
1789       xmax = 10.;
1790       break;
1791
1792    case 17:
1793       label = "N-Sigma of pions in TOF";
1794       nbins = 200;
1795       xmin = -10.;
1796       xmax = 10.;
1797       break;
1798
1799    case 18:
1800       label = "N-Sigma of protons in TOF";
1801       nbins = 200;
1802       xmin = -10.;
1803       xmax = 10.;
1804       break;
1805
1806    case 19:
1807       label = "N-Sigma of kaons in TOF";
1808       nbins = 200;
1809       xmin = -10.;
1810       xmax = 10.;
1811       break;
1812
1813    case 20:
1814       label = "TPC PID determination";
1815       nbins = 5;
1816       xmin = 0.;
1817       xmax = 5.;
1818       break;
1819
1820    case 21:
1821       label = "ITS PID determination";
1822       nbins = 5;
1823       xmin = 0.;
1824       xmax = 5.;
1825       break;
1826
1827    case 22:
1828       label = "TOF PID determination";
1829       nbins = 5;
1830       xmin = 0.;
1831       xmax = 5.;
1832       break;
1833
1834    } // end of switch
1835 } // end of get dimension parameters PID
1836
1837 void AliAnalysisTaskEmcalJetHadEPpid::Terminate(Option_t *) {
1838   cout<<"#########################"<<endl;
1839   cout<<"#### DONE RUNNING!!! ####"<<endl;
1840   cout<<"#########################"<<endl;
1841 } // end of terminate
1842
1843 //________________________________________________________________________
1844 Int_t AliAnalysisTaskEmcalJetHadEPpid::AcceptMyJet(AliEmcalJet *jet) {
1845   //applies all jet cuts except pt
1846   if ((jet->Phi()<fPhimin)||(jet->Phi()>fPhimax)) return 0;
1847   if ((jet->Eta()<fEtamin)||(jet->Eta()>fEtamax)) return 0;
1848   if (jet->Area()<fAreacut) return 0;
1849   // prevents 0 area jets from sneaking by when area cut == 0
1850   if (jet->Area()==0) return 0;
1851   //exclude jets with extremely high pt tracks which are likely misreconstructed
1852   if(jet->MaxTrackPt()>100) return 0;
1853   //prevents 0 area jets from sneaking by when area cut == 0
1854   if (jet->Area()==0) return 0;
1855   //exclude jets with extremely high pt tracks which are likely misreconstructed
1856   if(jet->MaxTrackPt()>100) return 0;
1857
1858   //passed all above cuts
1859   return 1;
1860 }
1861
1862 //void AliAnalysisTaskEmcalJetHadEPpid::FillAnalysisSummaryHistogram() const
1863 void AliAnalysisTaskEmcalJetHadEPpid::SetfHistPIDcounterLabels(TH1* h) const
1864 {
1865     // fill the analysis summary histrogram, saves all relevant analysis settigns
1866     h->GetXaxis()->SetBinLabel(1, "TPC: Unidentified"); // 0.5
1867     h->GetXaxis()->SetBinLabel(2, "TPC: Pion"); // 1.5
1868     h->GetXaxis()->SetBinLabel(3, "TPC: Kaon"); // 2.5
1869     h->GetXaxis()->SetBinLabel(4, "TPC: Proton"); // 3.5
1870     h->GetXaxis()->SetBinLabel(5, "ITS: Unidentified"); // 4.5
1871     h->GetXaxis()->SetBinLabel(6, "ITS: Pion"); // 5.5
1872     h->GetXaxis()->SetBinLabel(7, "ITS: Kaon"); // 6.5
1873     h->GetXaxis()->SetBinLabel(8, "ITS: Proton"); // 7.5
1874     h->GetXaxis()->SetBinLabel(9, "TOF: Unidentified"); // 8.5
1875     h->GetXaxis()->SetBinLabel(10, "TOF: Pion"); // 9.5
1876     h->GetXaxis()->SetBinLabel(11, "TOF: Kaon"); // 10.5
1877     h->GetXaxis()->SetBinLabel(12, "TOF: Proton"); // 11.5
1878     h->GetXaxis()->SetBinLabel(14, "Unidentified tracks"); //13.5
1879
1880 }
1881
1882 //______________________________________________________________________
1883 THnSparse* AliAnalysisTaskEmcalJetHadEPpid::NewTHnSparseDCorr(const char* name, UInt_t entries) {
1884   // generate new THnSparseD, axes are defined in GetDimParamsD()
1885   Int_t count = 0;
1886   UInt_t tmp = entries;
1887   while(tmp!=0){
1888     count++;
1889     tmp = tmp &~ -tmp;  // clear lowest bit
1890   }
1891
1892   TString hnTitle(name);
1893   const Int_t dim = count;
1894   Int_t nbins[dim];
1895   Double_t xmin[dim];
1896   Double_t xmax[dim];
1897
1898   Int_t i=0;
1899   Int_t c=0;
1900   while(c<dim && i<32){
1901     if(entries&(1<<i)){
1902
1903       TString label("");
1904       GetDimParamsCorr(i, label, nbins[c], xmin[c], xmax[c]);
1905       hnTitle += Form(";%s",label.Data());
1906       c++;
1907     }
1908
1909     i++;
1910   }
1911   hnTitle += ";";
1912
1913   return new THnSparseD(name, hnTitle.Data(), dim, nbins, xmin, xmax);
1914 } // end of NewTHnSparseD
1915
1916 //______________________________________________________________________________________________
1917 void AliAnalysisTaskEmcalJetHadEPpid::GetDimParamsCorr(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
1918 {
1919    //stores label and binning of axis for THnSparse
1920    const Double_t pi = TMath::Pi();
1921
1922    switch(iEntry){
1923
1924     case 0:
1925       label = "V0 centrality (%)";
1926       nbins = 10;
1927       xmin = 0.;
1928       xmax = 100.;
1929       break;
1930
1931    case 1:
1932       label = "Jet p_{T}";
1933       nbins = 216;
1934       xmin = 0.;
1935       xmax = 216.;
1936       break;
1937
1938    case 2:
1939       label = "Relative angle: Jet and Reaction Plane";
1940       nbins = 48;
1941       xmin = 0.;
1942       xmax = 0.5*pi;
1943       break;
1944
1945    }// end of switch
1946 } // end of Correction (ME) sparse
1947
1948