PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetHadEPpid.cxx

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