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