PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetHadEPpid.cxx

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