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