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