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[] = {-50,-45,-40,-35,-30,-25,-20,-18,-16,-14,-12,-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10,12,14,16,18,20,25,30,35,40,45,50,60,70,80,90,100,120,140,160,180,200,250,300,350,400};
 578   //Double_t xlowtrPT[] = {0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.25,2.50,2.75,3.0,3.25,3.5,3.75,4.0,4.25,4.50,4.75,5.0,5.5,6.0,6.5,7.0,7.5,8.0,8.5,9.0,9.5,10.0,11.0,12.0,13.0,14.0,15.0,16.0,17.0,18.0,19.0,20.0,22.0,24.0,26.0,28.0,30.0,35.0,40.0,45.0,50.0,60.0,70.0,80.0,90.0,100.0};
 579   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};
 580   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};
 581
 582   // tracks: 51, jets: 26
 583   // number of bins you tell histogram should be (# in array - 1) because the last bin
 584   // is the right-most edge of the histogram
 585   // i.e. this is for PT and there are 57 numbers (bins) thus we have 56 bins in our histo
 586   Int_t nbinsjetPT = sizeof(xlowjetPT)/sizeof(Double_t) - 1;
 587   Int_t nbinstrPT = sizeof(xlowtrPT)/sizeof(Double_t) - 1;
 588
 589   // set up jet-hadron sparse
 590   UInt_t bitcoded = 0; // bit coded, see GetDimParams() below
 591   UInt_t cifras = 0;
 592   UInt_t bitcode = 0;  // bit coded, see GetDimParamsPID() below
 593   UInt_t bitcodeCorr = 0; // bit coded, see GetDimparamsCorr() below
 594   bitcoded = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<7 | 1<<8; // | 1<<9;
 595   //if(fDoEventMixing) {
 596     fhnJH = NewTHnSparseF("fhnJH", bitcoded);
 597
 598     if(dovarbinTHnSparse){
 599       fhnJH->GetAxis(1)->Set(nbinsjetPT, xlowjetPT);
 600       fhnJH->GetAxis(2)->Set(nbinstrPT, xlowtrPT);
 601     }
 602
 603     fOutput->Add(fhnJH);
 604   //}
 605
 606   bitcodeCorr = 1<<0 | 1<<1 | 1<<2 | 1<<3; // | 1<<4 | 1<<5;
 607   fhnCorr = NewTHnSparseFCorr("fhnCorr", bitcodeCorr);
 608   if(dovarbinTHnSparse) fhnCorr->GetAxis(1)->Set(nbinsjetPT, xlowjetPT);
 609   fOutput->Add(fhnCorr);
 610
 611   /*
 612     Double_t centralityBins[nCentralityBins+1];
 613     for(Int_t ic=0; ic<nCentralityBins+1; ic++){
 614       if(ic==nCentralityBins) centralityBins[ic]=500;
 615       else centralityBins[ic]=10.0*ic;
 616     }
 617   */
 618
 619   // set up centrality bins for mixed events
 620   // for pp we need mult bins for event mixing. Create binning here, to also make a histogram from it
 621   Int_t nCentralityBinspp = 8;
 622   //Double_t centralityBinspp[nCentralityBinspp+1];
 623   Double_t centralityBinspp[9] = {0.0, 4., 9, 15, 25, 35, 55, 100.0, 500.0};
 624
 625   // Setup for Pb-Pb collisions
 626   Int_t nCentralityBinsPbPb = 10; //100;
 627   Double_t centralityBinsPbPb[nCentralityBinsPbPb+1];
 628   for(Int_t ic=0; ic<nCentralityBinsPbPb; ic++){
 629       centralityBinsPbPb[ic]=10.0*ic; //1.0*ic;
 630   }
 631
 632   if(fBeam == 0) fHistMult = new TH1F("fHistMult","multiplicity",nCentralityBinspp,centralityBinspp);
 633   if(fBeam == 1) fHistMult = new TH1F("fHistMult","multiplicity",nCentralityBinsPbPb,centralityBinsPbPb);
 634   //if(AliAnalysisTaskEmcal::GetBeamType() == 0) fHistMult = new TH1F("fHistMult","multiplicity",nCentralityBinspp,centralityBinspp);
 635   //if(AliAnalysisTaskEmcal::GetBeamType() == 1) fHistMult = new TH1F("fHistMult","multiplicity",nCentralityBinsPbPb,centralityBinsPbPb);
 636 //  fOutput->Add(fHistMult);
 637
 638   // Event Mixing
 639   Int_t trackDepth = fMixingTracks;
 640   Int_t poolsize   = 1000;  // Maximum number of events, ignored in the present implemented of AliEventPoolManager
 641   Int_t nZvtxBins  = 5+1+5;
 642   Double_t vertexBins[] = { -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10};
 643   Double_t* zvtxbin = vertexBins;
 644   if(fBeam == 0) fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBinspp, centralityBinspp, nZvtxBins, zvtxbin);
 645   if(fBeam == 1) fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBinsPbPb, centralityBinsPbPb, nZvtxBins, zvtxbin);
 646 //  if(GetBeamType() == 0) fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBinspp, centralityBinspp, nZvtxBins, zvtxbin);
 647 //  if(GetBeamType() == 1) fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBinsPbPb, centralityBinsPbPb, nZvtxBins, zvtxbin);
 648
 649   // set up event mixing sparse
 650   if(fDoEventMixing){
 651     cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<7 | 1<<8; // | 1<<9;
 652     fhnMixedEvents = NewTHnSparseF("fhnMixedEvents", cifras);
 653
 654     if(dovarbinTHnSparse){
 655      fhnMixedEvents->GetAxis(1)->Set(nbinsjetPT, xlowjetPT);
 656      fhnMixedEvents->GetAxis(2)->Set(nbinstrPT, xlowtrPT);
 657     }
 658
 659     fOutput->Add(fhnMixedEvents);
 660   } // end of do-eventmixing
 661
 662   // set up PID sparse
 663   if(doPID){
 664     // ****************************** PID *****************************************************
 665     // set up PID handler
 666     AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 667     AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
 668     if(!inputHandler) {
 669         AliFatal("Input handler needed");
 670                 return;
 671     }
 672
 673     // PID response object
 674     fPIDResponse = inputHandler->GetPIDResponse();
 675     if (!fPIDResponse) {
 676         AliError("PIDResponse object was not created");
 677             return;
 678     }
 679     // *****************************************************************************************
 680
 681     // PID counter
 682     fHistPID = new TH1F("fHistPID","PID Counter", 15, 0.5, 15.5);
 683     SetfHistPIDcounterLabels(fHistPID);
 684     fOutput->Add(fHistPID);
 685
 686     if(allpidAXIS) {
 687       bitcode = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<7 | 1<<8 | 1<<9 |
 688               1<<10 | 1<<11 | 1<<12 | 1<<13 | 1<<14 | 1<<15 | 1<<16 | 1<<17 | 1<<18 | 1<<19 |
 689               1<<20;
 690       fhnPID = NewTHnSparseFPID("fhnPID", bitcode);
 691     } else {
 692           bitcode = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<7 | 1<<8 | 1<<9 |
 693               1<<10 | 1<<11 | 1<<12 | 1<<13;
 694       fhnPID = NewTHnSparseFPID("fhnPID", bitcode);
 695         }
 696
 697     if(dovarbinTHnSparse){
 698      fhnPID->GetAxis(1)->Set(nbinstrPT, xlowtrPT);
 699      fhnPID->GetAxis(8)->Set(nbinsjetPT, xlowjetPT);
 700     }
 701
 702     fOutput->Add(fhnPID);
 703   } // end of do-PID
 704
 705   // =========== Switch on Sumw2 for all histos ===========
 706   for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
 707     TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
 708     if (h1){
 709       h1->Sumw2();
 710       continue;
 711     }
 712     TH2 *h2 = dynamic_cast<TH2*>(fOutput->At(i));
 713     if (h2){
 714       h2->Sumw2();
 715       continue;
 716     }
 717     TH3 *h3 = dynamic_cast<TH3*>(fOutput->At(i));
 718     if (h3){
 719       h3->Sumw2();
 720       continue;
 721     }
 722     THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
 723     if(hn)hn->Sumw2();
 724   }
 725
 726   PostData(1, fOutput);
 727 }
 728
 729 //_________________________________________________________________________
 730 void AliAnalysisTaskEmcalJetHadEPpid::ExecOnce()
 731 {
 732   AliAnalysisTaskEmcalJet::ExecOnce();
 733
 734   if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
 735   if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
 736   if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
 737 }
 738
 739 //_________________________________________________________________________
 740 Bool_t AliAnalysisTaskEmcalJetHadEPpid::Run()
 741 { // Main loop called for each event
 742   // TEST TEST TEST TEST for OBJECTS!
 743
 744   fHistEventQA->Fill(1); // All Events that get entered
 745
 746   UInt_t trig = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
 747   if(trig == 0) fHistEventSelectionQA->Fill(1);
 748   if(trig & AliVEvent::kAny) fHistEventSelectionQA->Fill(2);
 749   if(trig & AliVEvent::kAnyINT) fHistEventSelectionQA->Fill(3);
 750   if(trig & AliVEvent::kMB) fHistEventSelectionQA->Fill(4);
 751   if(trig & AliVEvent::kINT7) fHistEventSelectionQA->Fill(5);
 752   if(trig & AliVEvent::kEMC1) fHistEventSelectionQA->Fill(6);
 753   if(trig & AliVEvent::kEMC7) fHistEventSelectionQA->Fill(7);
 754   if(trig & AliVEvent::kEMC8) fHistEventSelectionQA->Fill(8);
 755   if(trig & AliVEvent::kEMCEJE) fHistEventSelectionQA->Fill(9);
 756   if(trig & AliVEvent::kEMCEGA) fHistEventSelectionQA->Fill(10);
 757   if(trig & AliVEvent::kCentral) fHistEventSelectionQA->Fill(11);
 758   if(trig & AliVEvent::kSemiCentral) fHistEventSelectionQA->Fill(12);
 759   if(trig & AliVEvent::kINT8) fHistEventSelectionQA->Fill(13);
 760
 761   if(trig & (AliVEvent::kEMCEJE | AliVEvent::kMB)) fHistEventSelectionQA->Fill(14);
 762   if(trig & (AliVEvent::kEMCEGA | AliVEvent::kMB)) fHistEventSelectionQA->Fill(15);
 763   if(trig & (AliVEvent::kAnyINT | AliVEvent::kMB)) fHistEventSelectionQA->Fill(16);
 764
 765   if(trig & (AliVEvent::kEMCEJE & AliVEvent::kMB)) fHistEventSelectionQA->Fill(17);
 766   if(trig & (AliVEvent::kEMCEGA & AliVEvent::kMB)) fHistEventSelectionQA->Fill(18);
 767   if(trig & (AliVEvent::kAnyINT & AliVEvent::kMB)) fHistEventSelectionQA->Fill(19);
 768
 769   if(GetBeamType() == 1) {
 770     if(!fLocalRho){
 771       AliError(Form("Couldn't get fLocalRho object, try to get it from Event based on name\n"));
 772       fLocalRho = GetLocalRhoFromEvent(fLocalRhoName);
 773       if(!fLocalRho) return kTRUE;
 774     }
 775   } // check for LocalRho object if PbPb data
 776
 777   if(!fTracks){
 778     AliError(Form("No fTracks object!!\n"));
 779     return kTRUE;
 780   }
 781   if(!fJets){
 782     AliError(Form("No fJets object!!\n"));
 783     return kTRUE;
 784   }
 785
 786   fHistEventQA->Fill(2); // events after object check
 787
 788   // what kind of event do we have: AOD or ESD?
 789   Bool_t useAOD;
 790   if (dynamic_cast<AliAODEvent*>(InputEvent())) useAOD = kTRUE;
 791   else useAOD = kFALSE;
 792
 793   // if we have ESD event, set up ESD object
 794   if(!useAOD){
 795     fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 796     if (!fESD) {
 797       AliError(Form("ERROR: fESD not available\n"));
 798       return kTRUE;
 799     }
 800   }
 801
 802   // if we have AOD event, set up AOD object
 803   if(useAOD){
 804     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 805     if(!fAOD) {
 806       AliError(Form("ERROR: fAOD not available\n"));
 807       return kTRUE;
 808     }
 809   }
 810
 811   fHistEventQA->Fill(3); // events after Aod/esd check
 812
 813   // get centrality
 814   Int_t centbin = GetCentBin(fCent);
 815   if (makeQAhistos) fHistCentrality->Fill(fCent); // won't be filled in pp collision (Keep this in mind!)
 816
 817   // BEAM TYPE enumerator: kNA = -1, kpp = 0, kAA = 1, kpA = 2
 818   // for pp analyses we will just use the first centrality bin
 819   if(GetBeamType() == 0) if (centbin == -1) centbin = 0;
 820   if(GetBeamType() == 1) if (centbin == -1) return kTRUE;
 821
 822   // if we are on PbPb data do cut on centrality > 90%, else by default DON'T
 823   if (GetBeamType() == 1) {
 824     // apply cut to event on Centrality > 90%
 825     if(fCent>90) return kTRUE;
 826   }
 827
 828   fHistEventQA->Fill(4);  // events after centrality check
 829
 830   // get vertex information
 831   Double_t fvertex[3]={0,0,0};
 832   InputEvent()->GetPrimaryVertex()->GetXYZ(fvertex);
 833   Double_t zVtx=fvertex[2];
 834   if (makeQAhistos) fHistZvtx->Fill(zVtx);
 835
 836   // get z-vertex bin
 837   //Int_t zVbin = GetzVertexBin(zVtx);
 838
 839   // apply zVtx cut
 840   if(fabs(zVtx)>10.0) return kTRUE;
 841
 842   fHistEventQA->Fill(5); // events after zvertex check
 843
 844   // create pointer to list of input event
 845   TList *list = InputEvent()->GetList();
 846   if(!list) {
 847     AliError(Form("ERROR: list not attached\n"));
 848     return kTRUE;
 849   }
 850
 851   fHistEventQA->Fill(6); // events after list check
 852
 853   // initialize TClonesArray pointers to jets and tracks
 854   TClonesArray *jets = 0;
 855   TClonesArray *tracks = 0;
 856   TClonesArray *tracksME = 0;
 857
 858   // get Tracks object
 859   tracks = dynamic_cast<TClonesArray*>(list->FindObject(fTracks));
 860   if (!tracks) {
 861     AliError(Form("Pointer to tracks %s == 0", fTracks->GetName()));
 862     return kTRUE;
 863   } // verify existence of tracks
 864
 865   // get ME Tracks object
 866   tracksME = dynamic_cast<TClonesArray*>(list->FindObject(fTracksNameME));
 867   if (!tracksME) {
 868     AliError(Form("Pointer to tracks %s == 0", fTracksNameME.Data()));
 869     return kTRUE;
 870   } // verify existence of tracks
 871
 872   // get Jets object
 873   jets = dynamic_cast<TClonesArray*>(list->FindObject(fJets));
 874   if(!jets){
 875     AliError(Form("Pointer to jets %s == 0", fJets->GetName()));
 876     return kTRUE;
 877   } // verify existence of jets
 878
 879   fHistEventQA->Fill(7);  // events after track/jet pointer check
 880
 881   // get number of jets and tracks
 882   const Int_t Njets = jets->GetEntries();
 883   const Int_t Ntracks = tracks->GetEntries();
 884   if(Ntracks<1)   return kTRUE;
 885   if(Njets<1)     return kTRUE;
 886
 887   fHistEventQA->Fill(8); // events after #track and jets < 1 check
 888
 889   if (makeQAhistos) fHistMult->Fill(Ntracks);  // fill multiplicity distribution
 890
 891   // initialize track parameters
 892   Int_t iTT=-1;
 893   Double_t ptmax=-10;
 894   Int_t NtrackAcc = 0;
 895
 896   fVevent = dynamic_cast<AliVEvent*>(InputEvent());
 897   if (!fVevent) {
 898     printf("ERROR: fVEvent not available\n");
 899     return kTRUE;
 900   }
 901
 902   // fill event mixing QA
 903   if(trig & AliVEvent::kEMCEGA) fHistCentZvertGA->Fill(fCent, zVtx);
 904   if(trig & AliVEvent::kEMCEJE) fHistCentZvertJE->Fill(fCent, zVtx);
 905   if(trig & AliVEvent::kMB) fHistCentZvertMB->Fill(fCent, zVtx);
 906   if(trig & AliVEvent::kAny) fHistCentZvertAny->Fill(fCent, zVtx);
 907
 908   // loop over tracks - to get hardest track (highest pt)
 909   for (Int_t iTracks = 0; iTracks < Ntracks; iTracks++){
 910     AliVParticle* Vtrack = static_cast<AliVParticle*>(tracks->At(iTracks));
 911     if (!Vtrack) {
 912       printf("ERROR: Could not receive track %d\n", iTracks);
 913       continue;
 914     }
 915
 916     AliVTrack *track = dynamic_cast<AliVTrack*>(Vtrack);
 917     if(!track) continue;
 918
 919     // track cuts
 920     if(TMath::Abs(track->Eta())>0.9) continue;
 921     if(track->Pt()<0.15) continue;
 922     //iCount++;
 923     NtrackAcc++;
 924
 925     if(track->Pt()>ptmax){
 926       ptmax=track->Pt();             // max pT track
 927       iTT=iTracks;                   // trigger tracks
 928     } // check if Pt>maxpt
 929
 930     if (makeQAhistos) fHistTrackPhi->Fill(track->Phi());
 931     if (makeQAhistos) fHistTrackPt[centbin]->Fill(track->Pt());
 932     if (makeQAhistos) fHistTrackPtallcent->Fill(track->Pt());
 933   } // end of loop over tracks
 934
 935   // get rho from event and fill relative histo's
 936   fRho = GetRhoFromEvent(fRhoName);
 937   fRhoVal = fRho->GetVal();
 938
 939   if (makeQAhistos) {
 940     fHistRhovsdEP[centbin]->Fill(fRhoVal,fEPV0); // Global Rho vs delta Event Plane angle
 941     fHistRhovsCent->Fill(fCent,fRhoVal);        // Global Rho vs Centrality
 942     fHistEP0[centbin]->Fill(fEPV0);
 943     fHistEP0A[centbin]->Fill(fEPV0A);
 944     fHistEP0C[centbin]->Fill(fEPV0C);
 945     fHistEPAvsC[centbin]->Fill(fEPV0A,fEPV0C);
 946   }
 947
 948   // initialize jet parameters
 949   Int_t ijethi=-1;
 950   Double_t highestjetpt=0.0;
 951   Int_t passedTTcut=0;
 952   Int_t NjetAcc = 0;
 953   Double_t leadhadronPT = 0;
 954
 955   // loop over jets in an event - to find highest jet pT and apply some cuts
 956   for (Int_t ijet = 0; ijet < Njets; ijet++){
 957     // get our jets
 958     AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(ijet));
 959     if (!jet) continue;
 960
 961     // apply jet cuts
 962     if ((jet->Phi()<fPhimin)||(jet->Phi()>fPhimax)) continue;
 963     if ((jet->Eta()<fEtamin)||(jet->Eta()>fEtamax)) continue;
 964     if (makeQAhistos) fHistAreavsRawPt[centbin]->Fill(jet->Pt(),jet->Area());
 965     if(!AcceptMyJet(jet)) continue;
 966
 967     NjetAcc++;                     // # of accepted jets
 968
 969     // if FlavourJetAnalysis, get desired FlavTag and check against Jet
 970     if(doFlavourJetAnalysis) { if(!AcceptFlavourJet(jet, fJetFlavTag)) continue;}
 971
 972     // use this to get total # of jets passing cuts in events!!!!!!!!
 973     if (makeQAhistos) fHistJetPhi->Fill(jet->Phi()); // Jet Phi histogram (filled)
 974
 975     // get highest Pt jet in event
 976     if(highestjetpt<jet->Pt()){
 977       ijethi=ijet;
 978       highestjetpt=jet->Pt();
 979     }
 980   } // end of looping over jets
 981
 982   // accepted jets
 983   fHistNjetvsCent->Fill(fCent,NjetAcc);
 984   Int_t NJETAcc = 0;
 985   fHistEventQA->Fill(9); // events after track/jet loop to get highest pt
 986
 987   // loop over jets in event and make appropriate cuts
 988   for (Int_t ijet = 0; ijet < Njets; ++ijet) {
 989      AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(ijet));
 990      if (!jet) continue;
 991
 992      if (!(trig & fTriggerEventType)) continue;
 993
 994          // (should probably be higher..., but makes a cut on jet pT)
 995      if (jet->Pt()<0.1) continue;
 996      // do we accept jet? apply jet cuts
 997      if (!AcceptMyJet(jet)) continue;
 998
 999      // if FlavourJetAnalysis, get desired FlavTag and check against Jet
1000      if(doFlavourJetAnalysis) { if(!AcceptFlavourJet(jet, fJetFlavTag)) continue;}
1001
1002      fHistEventQA->Fill(10); // accepted jets
1003
1004      // check on lead jet
1005      Double_t leadjet=0;
1006      if (ijet==ijethi) leadjet=1;
1007
1008      // check on leading hadron pt
1009      if (ijet==ijethi) leadhadronPT = GetLeadingHadronPt(jet);
1010
1011      // initialize and calculate various parameters: pt, eta, phi, rho, etc...
1012      Double_t jetphi = jet->Phi();      // phi of jet
1013      NJETAcc++;   // # accepted jets
1014      Double_t jeteta = jet->Eta();     // ETA of jet
1015      Double_t jetPt = -500;
1016      Double_t jetPtGlobal = -500;
1017      Double_t jetPtLocal = -500;            // initialize corr jet pt
1018      if(GetBeamType() == 1) {
1019        fLocalRhoVal = fLocalRho->GetLocalVal(jetphi, fJetRad); //GetJetRadius(0)); // get local rho value
1020        jetPtLocal = jet->Pt()-jet->Area()*fLocalRhoVal; // corrected pT of jet using Rho modulated for V2 and V3
1021      }
1022      jetPt = jet->Pt();
1023      jetPtGlobal = jet->Pt()-jet->Area()*fRhoVal;  // corrected pT of jet from rho value
1024      Double_t dEP = -500;                    // initialize angle between jet and event plane
1025      dEP = RelativeEPJET(jetphi,fEPV0); // angle betweeen jet and event plane
1026
1027      // make histo's
1028      if(makeQAhistos) fHistJetPtvsTrackPt[centbin]->Fill(jetPt,jet->MaxTrackPt());
1029      if(makeQAhistos) fHistRawJetPtvsTrackPt[centbin]->Fill(jetPt,jet->MaxTrackPt());
1030      if(makeQAhistos) fHistJetPtcorrGlRho[centbin]->Fill(jetPtGlobal);
1031      if(makeQAhistos) fHistJetPtvsdEP[centbin]->Fill(jetPt, dEP);
1032      if(makeQAhistos) fHistJetEtaPhiPt[centbin]->Fill(jetPt,jet->Eta(),jet->Phi());
1033      if(makeQAhistos) fHistJetPtArea[centbin]->Fill(jetPt,jet->Area());
1034      if(makeQAhistos) fHistJetEtaPhi->Fill(jet->Eta(),jet->Phi());     // fill jet eta-phi distribution histo
1035      if(makeextraCORRhistos) fHistJetPtNcon[centbin]->Fill(jetPt,jet->GetNumberOfConstituents());
1036      if(makeextraCORRhistos) fHistJetPtNconCh[centbin]->Fill(jetPt,jet->GetNumberOfTracks());
1037      if(makeextraCORRhistos) fHistJetPtNconEm[centbin]->Fill(jetPt,jet->GetNumberOfClusters());
1038      if (makeoldJEThadhistos) fHistJetPt[centbin]->Fill(jet->Pt());  // fill #jets vs pT histo
1039      //fHistDeltaPtvsArea->Fill(jetPt,jet->Area());
1040
1041      // make histo's with BIAS applied
1042      if (jet->MaxTrackPt()>fTrkBias){
1043        if(makeBIAShistos) fHistJetPtvsdEPBias[centbin]->Fill(jetPt, dEP);
1044        if(makeBIAShistos) fHistJetEtaPhiPtBias[centbin]->Fill(jetPt,jet->Eta(),jet->Phi());
1045        if(makeextraCORRhistos) fHistJetPtAreaBias[centbin]->Fill(jetPt,jet->Area());
1046        if(makeextraCORRhistos) fHistJetPtNconBias[centbin]->Fill(jetPt,jet->GetNumberOfConstituents());
1047        if(makeextraCORRhistos) fHistJetPtNconBiasCh[centbin]->Fill(jetPt,jet->GetNumberOfTracks());
1048        if(makeextraCORRhistos) fHistJetPtNconBiasEm[centbin]->Fill(jetPt,jet->GetNumberOfClusters());
1049      }
1050
1051     //if(leadjet && centbin==0){
1052     //  if(makeextraCORRhistos) fHistJetPt[centbin+1]->Fill(jet->Pt());
1053     //}
1054     if ((jet->MaxTrackPt()>fTrkBias) || (jet->MaxClusterPt()>fClusBias)){
1055       if (makeoldJEThadhistos){
1056         fHistJetPtBias[centbin]->Fill(jet->Pt());
1057         //if(leadjet && centbin==0) fHistJetPtBias[centbin+1]->Fill(jet->Pt());
1058       }
1059     }  // end of MaxTrackPt>ftrkBias or maxclusterPt>fclusBias
1060
1061     // do we have trigger tracks
1062     if(iTT>0){
1063       AliVTrack* TT = static_cast<AliVTrack*>(tracks->At(iTT));
1064       if(TMath::Abs(jet->Phi()-TT->Phi()-TMath::Pi())<0.6) passedTTcut=1;
1065       else passedTTcut=0;
1066     } // end of check on iTT > 0
1067
1068     if(passedTTcut) {
1069       if (makeoldJEThadhistos) fHistJetPtTT[centbin]->Fill(jet->Pt());
1070     }
1071
1072     // cut on HIGHEST jet pt in event (15 GeV default)
1073     //if (highestjetpt>fJetPtcut) {
1074     if (jet->Pt() > fJetPtcut) {
1075       fHistEventQA->Fill(11); // jets meeting pt threshold
1076
1077       // does our max track or cluster pass the bias?
1078       if ((jet->MaxTrackPt()>fTrkBias) || (jet->MaxClusterPt()>fClusBias)){
1079         // set up and fill Jet-Hadron Correction THnSparse
1080         Double_t CorrEntries[4] = {fCent, jet->Pt(), dEP, zVtx};
1081         fhnCorr->Fill(CorrEntries);    // fill Sparse Histo with Correction entries
1082
1083         if(GetBeamType() == 1) fHistLocalRhoJetpt->Fill(jetPtLocal);
1084       }
1085
1086       // loop over all track for an event containing a jet with a pT>fJetPtCut  (15)GeV
1087       for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
1088         AliVParticle* Vtrackass = static_cast<AliVParticle*>(tracks->At(iTracks));
1089         if (!Vtrackass) {
1090           printf("ERROR: Could not receive track %d\n", iTracks);
1091           continue;
1092         }
1093
1094         AliVTrack *track = dynamic_cast<AliVTrack*>(Vtrackass);
1095         if (!track) {
1096           AliError(Form("Couldn't get AliVtrack %d\n", iTracks));
1097           continue;
1098         }
1099
1100         // apply track cuts
1101         if(TMath::Abs(track->Eta())>fTrkEta) continue;
1102         if (track->Pt()<0.15) continue;
1103
1104         fHistEventQA->Fill(12); // accepted tracks in events from trigger jets
1105
1106         // calculate and get some track parameters
1107                 Double_t trCharge = -99;
1108         trCharge = track->Charge();
1109         Double_t tracketa=track->Eta();   // eta of track
1110         Double_t deta=tracketa-jeteta;    // dETA between track and jet
1111         //Double_t dR=sqrt(deta*deta+dphijh*dphijh);     // difference of R between jet and hadron track
1112
1113         Int_t ieta = -1;       // initialize deta bin
1114         Int_t iptjet = -1;     // initialize jet pT bin
1115         if (makeoldJEThadhistos)  {
1116                   ieta=GetEtaBin(deta);             // bin of eta
1117               if(ieta<0) continue;              // double check we don't have a negative array index
1118           iptjet=GetpTjetBin(jet->Pt());    // bin of jet pt
1119               if(iptjet<0) continue;                      // double check we don't have a negative array index
1120                 }
1121
1122         // dPHI between jet and hadron
1123         Double_t dphijh = RelativePhi(jet->Phi(), track->Phi()); // angle between jet and hadron
1124
1125         // fill some jet-hadron histo's
1126         if (makeoldJEThadhistos) fHistJetH[centbin][iptjet][ieta]->Fill(dphijh,track->Pt());  // fill jet-hadron dPHI--track pT distribution
1127         if(makeQAhistos) fHistJetHEtaPhi->Fill(deta,dphijh);                          // fill jet-hadron  eta--phi distribution
1128             fHistJetHaddPHI->Fill(dphijh);
1129         if(passedTTcut){
1130           if (makeoldJEThadhistos) fHistJetHTT[centbin][iptjet][ieta]->Fill(dphijh,track->Pt());
1131         }
1132
1133         // does our max track or cluster pass the bias?
1134         if ((jet->MaxTrackPt()>fTrkBias) || (jet->MaxClusterPt()>fClusBias)){
1135           // set up and fill Jet-Hadron THnSparse
1136           Double_t triggerEntries[9] = {fCent, jet->Pt(), track->Pt(), deta, dphijh, dEP, zVtx, trCharge, leadjet};
1137           fhnJH->Fill(triggerEntries);    // fill Sparse Histo with trigger entries
1138
1139           // fill histo's
1140           if(makeQAhistos) fHistSEphieta->Fill(dphijh, deta); // single event distribution
1141           if (makeoldJEThadhistos) fHistJetHBias[centbin][iptjet][ieta]->Fill(dphijh,track->Pt());
1142
1143           if (makeBIAShistos) {
1144             fHistJetHaddPhiBias->Fill(dphijh);
1145
1146             // in plane and out of plane histo's
1147             if( dEP>0 && dEP<=(TMath::Pi()/6) ){
1148               // we are IN plane
1149               fHistJetHaddPhiINBias->Fill(dphijh);
1150             }else if( dEP>(TMath::Pi()/3) && dEP<=(TMath::Pi()/2) ){
1151               // we are OUT of PLANE
1152               fHistJetHaddPhiOUTBias->Fill(dphijh);
1153             }else if( dEP>(TMath::Pi()/6) && dEP<=(TMath::Pi()/3) ){
1154               // we are in middle of plane
1155               fHistJetHaddPhiMIDBias->Fill(dphijh);
1156             }
1157                   } // BIAS histos switch
1158         } // end of check maxtrackpt>ftrackbias or maxclusterpt>fclustbias
1159
1160         // **************************************************************************************************************
1161         // *********************************** PID **********************************************************************
1162         // **************************************************************************************************************
1163         if(doPIDtrackBIAS){
1164           //if(ptmax < fTrkBias) continue;    // force PID to happen when max track pt > 5.0 GeV
1165           if(leadhadronPT < fTrkBias) continue; // force PID to happen when lead hadron pt > 5.0 GeV
1166         }
1167
1168         // some variables for PID
1169         Double_t pt = -999;
1170         Double_t dEdx = -999;
1171         Double_t ITSsig = -999;
1172         Double_t TOFsig = -999;
1173         Double_t charge = -999;
1174
1175         // nSigma of particles in TPC, TOF, and ITS
1176         Double_t nSigmaPion_TPC, nSigmaProton_TPC, nSigmaKaon_TPC;
1177         Double_t nSigmaPion_TOF, nSigmaProton_TOF, nSigmaKaon_TOF;
1178         Double_t nSigmaPion_ITS, nSigmaProton_ITS, nSigmaKaon_ITS;
1179
1180         if(doPID){
1181           // get parameters of track
1182           charge = track->Charge();    // charge of track
1183           pt     = track->Pt();        // pT of track
1184
1185           // extra attempt
1186           AliVEvent *vevent=InputEvent();
1187           if (!vevent||!fPIDResponse) return kTRUE; // just return, maybe put at beginning
1188
1189           fHistEventQA->Fill(13); // check for AliVEvent and fPIDresponse objects
1190
1191 ///////////////////////////////////////
1192
1193           // get detector signals
1194           dEdx = track->GetTPCsignal();
1195           ITSsig = track->GetITSsignal();
1196           TOFsig = track->GetTOFsignal();
1197
1198           // TPC nSigma's
1199           nSigmaPion_TPC = fPIDResponse->NumberOfSigmasTPC(track,AliPID::kPion);
1200           nSigmaKaon_TPC = fPIDResponse->NumberOfSigmasTPC(track,AliPID::kKaon);
1201           nSigmaProton_TPC = fPIDResponse->NumberOfSigmasTPC(track,AliPID::kProton);
1202
1203           // TOF nSigma's
1204           nSigmaPion_TOF = fPIDResponse->NumberOfSigmasTOF(track,AliPID::kPion);
1205           nSigmaKaon_TOF = fPIDResponse->NumberOfSigmasTOF(track,AliPID::kKaon);
1206           nSigmaProton_TOF = fPIDResponse->NumberOfSigmasTOF(track,AliPID::kProton);
1207
1208           // ITS nSigma's
1209           nSigmaPion_ITS = fPIDResponse->NumberOfSigmasITS(track,AliPID::kPion);
1210           nSigmaKaon_ITS = fPIDResponse->NumberOfSigmasITS(track,AliPID::kKaon);
1211           nSigmaProton_ITS = fPIDResponse->NumberOfSigmasITS(track,AliPID::kProton);
1212
1213           // fill detector signal histograms
1214           if (makeQAhistos) fHistTPCdEdX->Fill(pt, dEdx);
1215           if (makeQAhistos) fHistITSsignal->Fill(pt, ITSsig);
1216           //if (makeQAhistos) fHistTOFsignal->Fill(pt, TOFsig);
1217
1218           // Tests to PID pions, kaons, and protons,          (default is undentified tracks)
1219           Double_t nPIDtpc = 0;
1220           Double_t nPIDits = 0;
1221           Double_t nPIDtof = 0;
1222           Double_t nPID = -99;
1223
1224           // check track has pT < 0.900 GeV  - use TPC pid
1225           if (pt<0.900 && dEdx>0) {
1226                 nPIDtpc = 4;
1227             nPID = 1;
1228
1229             // PION check - TPC
1230             if (TMath::Abs(nSigmaPion_TPC)<2 && TMath::Abs(nSigmaKaon_TPC)>2 && TMath::Abs(nSigmaProton_TPC)>2 ){
1231               isPItpc = kTRUE;
1232               nPIDtpc = 1;
1233               nPID=2;
1234             }else isPItpc = kFALSE;
1235
1236             // KAON check - TPC
1237             if (TMath::Abs(nSigmaKaon_TPC)<2 && TMath::Abs(nSigmaPion_TPC)>3 && TMath::Abs(nSigmaProton_TPC)>2 ){
1238               isKtpc = kTRUE;
1239               nPIDtpc = 2;
1240               nPID=3;
1241             }else isKtpc = kFALSE;
1242
1243             // PROTON check - TPC
1244             if (TMath::Abs(nSigmaProton_TPC)<2 && TMath::Abs(nSigmaPion_TPC)>3 && TMath::Abs(nSigmaKaon_TPC)>2 ){
1245               isPtpc = kTRUE;
1246               nPIDtpc = 3;
1247               nPID=4;
1248             }else isPtpc = kFALSE;
1249           }  // cut on track pT for TPC
1250
1251           // check track has pT < 0.500 GeV - use ITS pid
1252           if (pt<0.500 && ITSsig>0) {
1253             nPIDits = 4;
1254             nPID = 5;
1255
1256             // PION check - ITS
1257             if (TMath::Abs(nSigmaPion_ITS)<2 && TMath::Abs(nSigmaKaon_ITS)>2 && TMath::Abs(nSigmaProton_ITS)>2 ){
1258               isPIits = kTRUE;
1259               nPIDits = 1;
1260                   nPID=6;
1261             }else isPIits = kFALSE;
1262
1263             // KAON check - ITS
1264             if (TMath::Abs(nSigmaKaon_ITS)<2 && TMath::Abs(nSigmaPion_ITS)>3 && TMath::Abs(nSigmaProton_ITS)>2 ){
1265               isKits = kTRUE;
1266               nPIDits = 2;
1267               nPID=7;
1268             }else isKits = kFALSE;
1269
1270             // PROTON check - ITS
1271             if (TMath::Abs(nSigmaProton_ITS)<2 && TMath::Abs(nSigmaPion_ITS)>3 && TMath::Abs(nSigmaKaon_ITS)>2 ){
1272               isPits = kTRUE;
1273               nPIDits = 3;
1274                   nPID=8;
1275             }else isPits = kFALSE;
1276           }  // cut on track pT for ITS
1277
1278           // check track has 0.900 GeV < pT < 2.500 GeV - use TOF pid
1279           if (pt>0.900 && pt<2.500 && TOFsig>0) {
1280                 nPIDtof = 4;
1281             nPID = 9;
1282
1283             // PION check - TOF
1284             if (TMath::Abs(nSigmaPion_TOF)<2 && TMath::Abs(nSigmaKaon_TOF)>2 && TMath::Abs(nSigmaProton_TOF)>2 ){
1285               isPItof = kTRUE;
1286               nPIDtof = 1;
1287               nPID=10;
1288             }else isPItof = kFALSE;
1289
1290             // KAON check - TOF
1291             if (TMath::Abs(nSigmaKaon_TOF)<2 && TMath::Abs(nSigmaPion_TOF)>3 && TMath::Abs(nSigmaProton_TOF)>2 ){
1292               isKtof = kTRUE;
1293               nPIDtof = 2;
1294               nPID=11;
1295             }else isKtof = kFALSE;
1296
1297             // PROTON check - TOF
1298             if (TMath::Abs(nSigmaProton_TOF)<2 && TMath::Abs(nSigmaPion_TOF)>3 && TMath::Abs(nSigmaKaon_TOF)>2 ){
1299               isPtof = kTRUE;
1300               nPIDtof = 3;
1301               nPID=12;
1302             }else isPtof = kFALSE;
1303           }  // cut on track pT for TOF
1304
1305           if (nPID == -99) nPID = 14;
1306               fHistPID->Fill(nPID);
1307
1308           // PID sparse getting filled
1309           if (allpidAXIS) { // FILL ALL axis
1310                         Double_t pid_EntriesALL[21] = {fCent,pt,charge,deta,dphijh,leadjet,zVtx,dEP,jetPt,
1311                                       nSigmaPion_TPC, nSigmaPion_TOF, // pion nSig values in TPC/TOF
1312                                                                   nPIDtpc, nPIDits, nPIDtof,       // PID label for each detector
1313                                                                           nSigmaProton_TPC, nSigmaKaon_TPC,  // nSig in TPC
1314                                       nSigmaPion_ITS, nSigmaProton_ITS, nSigmaKaon_ITS,  // nSig in ITS
1315                                       nSigmaProton_TOF, nSigmaKaon_TOF,  // nSig in TOF
1316                                       }; //array for PID sparse
1317                     fhnPID->Fill(pid_EntriesALL);
1318                   } else {
1319             // PID sparse getting filled
1320             Double_t pid_Entries[14] = {fCent,pt,charge,deta,dphijh,leadjet,zVtx,dEP,jetPt,
1321                                       nSigmaPion_TPC, nSigmaPion_TOF, // pion nSig values in TPC/TOF
1322                                                                   nPIDtpc, nPIDits, nPIDtof       // PID label for each detector
1323                                       }; //array for PID sparse
1324             fhnPID->Fill(pid_Entries);   // fill Sparse histo of PID tracks
1325                   } // minimal pid sparse filling
1326
1327         } // end of doPID check
1328
1329             // get track pt bin
1330         Int_t itrackpt = -500;              // initialize track pT bin
1331         itrackpt = GetpTtrackBin(track->Pt());
1332
1333             // all tracks: jet hadron correlations in hadron pt bins
1334         if(makeextraCORRhistos) fHistJetHadbindPhi[itrackpt]->Fill(dphijh);
1335
1336         // in plane and out of plane jet-hadron histo's
1337         if( dEP>0 && dEP<=(TMath::Pi()/6) ){
1338           // we are IN plane
1339           if(makeextraCORRhistos) fHistJetHaddPhiINcent[centbin]->Fill(dphijh);
1340           fHistJetHaddPhiIN->Fill(dphijh);
1341           if(makeextraCORRhistos) fHistJetHadbindPhiIN[itrackpt]->Fill(dphijh);
1342           //fHistJetHaddPhiPtcentbinIN[itrackpt][centbin]->Fill(dphijh);
1343         }else if( dEP>(TMath::Pi()/3) && dEP<=(TMath::Pi()/2) ){
1344           // we are OUT of PLANE
1345           if(makeextraCORRhistos) fHistJetHaddPhiOUTcent[centbin]->Fill(dphijh);
1346           fHistJetHaddPhiOUT->Fill(dphijh);
1347           if(makeextraCORRhistos) fHistJetHadbindPhiOUT[itrackpt]->Fill(dphijh);
1348           //fHistJetHaddPhiPtcentbinOUT[itrackpt][centbin]->Fill(dphijh);
1349         }else if( dEP>(TMath::Pi()/6) && dEP<=(TMath::Pi()/3) ){
1350           // we are in the middle of plane
1351           if(makeextraCORRhistos) fHistJetHaddPhiMIDcent[centbin]->Fill(dphijh);
1352           fHistJetHaddPhiMID->Fill(dphijh);
1353           if(makeextraCORRhistos) fHistJetHadbindPhiMID[itrackpt]->Fill(dphijh);
1354         }
1355       } // loop over tracks found in event with highest JET pT > 10.0 GeV (change)
1356     } // jet pt cut
1357   } // jet loop
1358
1359   fHistEventQA->Fill(14); // events right before event mixing
1360
1361 // ***************************************************************************************************************
1362 // ******************************** Event MIXING *****************************************************************
1363 // ***************************************************************************************************************
1364
1365   // initialize object array of cloned picotracks
1366   TObjArray* tracksClone = 0x0;
1367
1368   // PbPb collisions - create cloned picotracks
1369   //if(GetBeamType() == 1) tracksClone = CloneAndReduceTrackList(tracks); // TEST
1370
1371   //Prepare to do event mixing
1372   if(fDoEventMixing>0){
1373     // event mixing
1374
1375     // 1. First get an event pool corresponding in mult (cent) and
1376     //    zvertex to the current event. Once initialized, the pool
1377     //    should contain nMix (reduced) events. This routine does not
1378     //    pre-scan the chain. The first several events of every chain
1379     //    will be skipped until the needed pools are filled to the
1380     //    specified depth. If the pool categories are not too rare, this
1381     //    should not be a problem. If they are rare, you could lose
1382     //    statistics.
1383
1384     // 2. Collect the whole pool's content of tracks into one TObjArray
1385     //    (bgTracks), which is effectively a single background super-event.
1386
1387     // 3. The reduced and bgTracks arrays must both be passed into
1388     //    FillCorrelations(). Also nMix should be passed in, so a weight
1389     //    of 1./nMix can be applied.
1390
1391     // mix jets from triggered events with tracks from MB events
1392     // get the trigger bit, need to change trigger bits between different runs
1393     UInt_t trigger = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
1394     // if event was not selected (triggered) for any reseason (should never happen) then return
1395     if (trigger==0)  return kTRUE;
1396
1397     // initialize event pools
1398     AliEventPool* pool = 0x0;
1399     AliEventPool* poolpp = 0x0;
1400     Double_t Ntrks = -999;
1401
1402     // pp collisions - get event pool
1403     if(GetBeamType() == 0) {
1404       Ntrks=(Double_t)Ntracks*1.0;
1405       //cout<<"Test.. Ntrks: "<<fPoolMgr->GetEventPool(Ntrks);
1406       poolpp = fPoolMgr->GetEventPool(Ntrks, zVtx); // for pp
1407     }
1408
1409     // PbPb collisions - get event pool
1410     if(GetBeamType() == 1) pool = fPoolMgr->GetEventPool(fCent, zVtx); // for PbPb? fcent
1411
1412     // if we don't have a pool, return
1413     if (!pool && !poolpp){
1414       if(GetBeamType() == 1) AliFatal(Form("No pool found for centrality = %f, zVtx = %f", fCent, zVtx));
1415       if(GetBeamType() == 0) AliFatal(Form("No pool found for multiplicity = %f, zVtx = %f", Ntrks, zVtx));
1416       return kTRUE;
1417     }
1418
1419     fHistEventQA->Fill(15); // mixed events cases that have pool
1420
1421     // initialize background tracks array
1422     TObjArray* bgTracks;
1423
1424     // next line might not apply for PbPb collisions
1425     // use only jets from EMCal-triggered events (for lhc11a use AliVEvent::kEMC1)
1426     //check for a trigger jet
1427     // fmixingtrack/10 ??
1428   if(GetBeamType() == 1) if(trigger & fTriggerEventType) { //kEMCEJE)) {
1429     if (pool->IsReady() || pool->NTracksInPool() > fNMIXtracks || pool->GetCurrentNEvents() >= fNMIXevents) {
1430
1431       // loop over jets (passing cuts?)
1432       for (Int_t ijet = 0; ijet < Njets; ijet++) {
1433         Double_t leadjet=0;
1434         if (ijet==ijethi) leadjet=1;
1435
1436         // get jet object
1437         AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(ijet));
1438             if (!jet) continue;
1439
1440                 // (should probably be higher..., but makes a cut on jet pT)
1441         if (jet->Pt()<0.1) continue;
1442         if (!AcceptMyJet(jet)) continue;
1443
1444         fHistEventQA->Fill(16); // event mixing jets
1445
1446         // set cut to do event mixing only if we have a jet meeting our pt threshold (bias applied below)
1447         if (jet->Pt()<fJetPtcut) continue;
1448
1449         // get number of current events in pool
1450         Int_t nMix = pool->GetCurrentNEvents();  // how many particles in pool to mix
1451
1452           // Fill for biased jet triggers only
1453         if ((jet->MaxTrackPt()>fTrkBias) || (jet->MaxClusterPt()>fClusBias)) {  // && jet->Pt() > fJetPtcut) {
1454           // Fill mixed-event histos here
1455           for (Int_t jMix=0; jMix<nMix; jMix++) {
1456                         fHistEventQA->Fill(17); // event mixing nMix
1457
1458                     // get jMix'th event
1459                         bgTracks = pool->GetEvent(jMix);
1460             const Int_t Nbgtrks = bgTracks->GetEntries();
1461             for(Int_t ibg=0; ibg<Nbgtrks; ibg++) {
1462               AliPicoTrack *part = static_cast<AliPicoTrack*>(bgTracks->At(ibg));
1463               if(!part) continue;
1464               if(TMath::Abs(part->Eta())>0.9) continue;
1465               if(part->Pt()<0.15) continue;
1466
1467               Double_t DEta = part->Eta()-jet->Eta();                // difference in eta
1468               Double_t DPhi = RelativePhi(jet->Phi(),part->Phi());   // difference in phi
1469               Double_t dEP = RelativeEPJET(jet->Phi(),fEPV0);        // difference between jet and EP
1470                       Double_t mixcharge = part->Charge();
1471               //Double_t DR=TMath::Sqrt(DPhi*DPhi+DEta*DEta);      // difference in R
1472
1473               // create / fill mixed event sparse
1474               Double_t triggerEntries[10] = {fCent,jet->Pt(),part->Pt(),DEta,DPhi,dEP,zVtx, mixcharge, leadjet}; //array for ME sparse
1475               fhnMixedEvents->Fill(triggerEntries,1./nMix);   // fill Sparse histo of mixed events
1476
1477               fHistEventQA->Fill(18); // event mixing - nbgtracks
1478               if(makeextraCORRhistos) fHistMEphieta->Fill(DPhi,DEta, 1./nMix);
1479             } // end of background track loop
1480           } // end of filling mixed-event histo's
1481         } // end of check for biased jet triggers
1482       } // end of jet loop
1483     } // end of check for pool being ready
1484   } // end EMC triggered loop
1485
1486 //=============================================================================================================
1487
1488     // use only jets from EMCal-triggered events (for lhc11a use AliVEvent::kEMC1)
1489 ///    if (trigger & AliVEvent::kEMC1) {
1490     // pp collisions
1491     if(GetBeamType() == 0) if(trigger & fTriggerEventType) { //kEMC1)) {
1492       if (poolpp->IsReady() || poolpp->NTracksInPool() > fNMIXtracks || poolpp->GetCurrentNEvents() >= fNMIXevents) {
1493
1494         // loop over jets (passing cuts?)
1495         for (Int_t ijet = 0; ijet < Njets; ijet++) {
1496           Double_t leadjet=0;
1497           if (ijet==ijethi) leadjet=1;
1498
1499           // get jet object
1500           AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(ijet));
1501                   if (!jet) continue;
1502
1503                   // (should probably be higher..., but makes a cut on jet pT)
1504           if (jet->Pt()<0.1) continue;
1505           if (!AcceptMyJet(jet)) continue;
1506
1507           fHistEventQA->Fill(16); // event mixing jets
1508
1509           // set cut to do event mixing only if we have a jet meeting our pt threshold (bias applied below)
1510                   if (jet->Pt()<fJetPtcut) continue;
1511
1512           // get number of current events in pool
1513           Int_t nMix = poolpp->GetCurrentNEvents();  // how many particles in pool to mix
1514
1515           // Fill for biased jet triggers only
1516           if ((jet->MaxTrackPt()>fTrkBias) || (jet->MaxClusterPt()>fClusBias)) {  // && jet->Pt() > fJetPtcut) {
1517             // Fill mixed-event histos here
1518             for (Int_t jMix=0; jMix<nMix; jMix++) {
1519                           fHistEventQA->Fill(17); // event mixing nMix
1520
1521                       // get jMix'th event
1522                           bgTracks = poolpp->GetEvent(jMix);
1523               const Int_t Nbgtrks = bgTracks->GetEntries();
1524               for(Int_t ibg=0; ibg<Nbgtrks; ibg++) {
1525                 AliPicoTrack *part = static_cast<AliPicoTrack*>(bgTracks->At(ibg));
1526                 if(!part) continue;
1527                 if(TMath::Abs(part->Eta())>0.9) continue;
1528                 if(part->Pt()<0.15) continue;
1529
1530                 Double_t DEta = part->Eta()-jet->Eta();                // difference in eta
1531                 Double_t DPhi = RelativePhi(jet->Phi(),part->Phi());   // difference in phi
1532                 Double_t dEP = RelativeEPJET(jet->Phi(),fEPV0);      // difference between jet and EP
1533                 Double_t mixcharge = part->Charge();
1534                 //Double_t DR=TMath::Sqrt(DPhi*DPhi+DEta*DEta);      // difference in R
1535
1536                 // create / fill mixed event sparse
1537                 Double_t triggerEntries[10] = {fCent,jet->Pt(),part->Pt(),DEta,DPhi,dEP,zVtx, mixcharge, leadjet}; //array for ME sparse
1538                 fhnMixedEvents->Fill(triggerEntries,1./nMix);   // fill Sparse histo of mixed events
1539
1540                 fHistEventQA->Fill(18); // event mixing - nbgtracks
1541                 if(makeextraCORRhistos) fHistMEphieta->Fill(DPhi,DEta, 1./nMix);
1542               } // end of background track loop
1543             } // end of filling mixed-event histo's
1544           } // end of check for biased jet triggers
1545         } // end of jet loop
1546       } // end of check for pool being ready
1547     } //end EMC triggered loop
1548
1549     // pp collisions
1550     // use only tracks from MB events (for lhc11a use AliVEvent::kMB)
1551 ///    if (trigger & AliVEvent::kMB) {
1552 //T    if (trigger & AliVEvent::kAnyINT){ // test
1553     if(GetBeamType() == 0) {
1554
1555       // use only tracks from MB events (for lhc11a use AliVEvent::kMB)
1556       if(trigger & fMixingEventType) { //kMB) {
1557
1558         // create a list of reduced objects. This speeds up processing and reduces memory consumption for the event pool
1559         tracksClone = CloneAndReduceTrackList(tracks);
1560
1561         // update pool if jet in event or not
1562         poolpp->UpdatePool(tracksClone);
1563
1564       } // check on track from MB events
1565     }
1566
1567     // PbPb collisions
1568     if(GetBeamType() == 1) {
1569
1570       // use only tracks from MB events
1571       if(trigger & fMixingEventType) { //kMB) {
1572
1573         // create a list of reduced objects. This speeds up processing and reduces memory consumption for the event pool
1574         tracksClone = CloneAndReduceTrackList(tracks);
1575
1576         // update pool if jet in event or not
1577         pool->UpdatePool(tracksClone);
1578
1579       } // MB events
1580     } // PbPb collisions
1581   } // end of event mixing
1582
1583   // print some stats on the event
1584   event++;
1585   fHistEventQA->Fill(19);  // events making it to end
1586
1587   if (doComments) {
1588     cout<<"Event #: "<<event<<"     Jet Radius: "<<fJetRad<<"     Constituent Pt Cut: "<<fConstituentCut<<endl;
1589     cout<<"# of jets: "<<Njets<<"      NjetAcc: "<<NjetAcc<<"      Highest jet pt: "<<highestjetpt<<"     leading hadron pt: "<<leadhadronPT<<endl;
1590     cout<<"# tracks: "<<Ntracks<<"      NtrackAcc: "<<NtrackAcc<<"      Highest track pt: "<<ptmax<<endl;
1591     cout<<" =============================================== "<<endl;
1592   }
1593
1594   return kTRUE;  // used when the function is of type bool
1595 }  // end of RUN
1596
1597 //________________________________________________________________________
1598 Int_t AliAnalysisTaskEmcalJetHadEPpid::GetCentBin(Double_t cent) const
1599 {  // Get centrality bin.
1600   Int_t centbin = -1;
1601   if (cent>=0 && cent<10)       centbin = 0;
1602   else if (cent>=10 && cent<20) centbin = 1;
1603   else if (cent>=20 && cent<30) centbin = 2;
1604   else if (cent>=30 && cent<40) centbin = 3;
1605   else if (cent>=40 && cent<50) centbin = 4;
1606   else if (cent>=50 && cent<90) centbin = 5;
1607
1608   return centbin;
1609 }
1610
1611 //________________________________________________________________________
1612 Double_t AliAnalysisTaskEmcalJetHadEPpid::RelativePhi(Double_t mphi,Double_t vphi) const
1613 { // function to calculate relative PHI
1614   double dphi = mphi-vphi;
1615
1616   // set dphi to operate on adjusted scale
1617   if(dphi<-0.5*TMath::Pi()) dphi+=2.*TMath::Pi();
1618   if(dphi>3./2.*TMath::Pi()) dphi-=2.*TMath::Pi();
1619
1620   // test
1621   if( dphi < -1.*TMath::Pi()/2 || dphi > 3.*TMath::Pi()/2 )
1622     AliWarning(Form("%s: dPHI not in range [-0.5*Pi, 1.5*Pi]!", GetName()));
1623
1624   return dphi; // dphi in [-0.5Pi, 1.5Pi]
1625 }
1626
1627 //_________________________________________________________________________
1628 Double_t AliAnalysisTaskEmcalJetHadEPpid:: RelativeEPJET(Double_t jetAng, Double_t EPAng) const
1629 { // function to calculate angle between jet and EP in the 1st quadrant (0,Pi/2)
1630   Double_t dphi = (EPAng - jetAng);
1631
1632   // ran into trouble with a few dEP<-Pi so trying this...
1633   if( dphi<-1*TMath::Pi() ){
1634     dphi = dphi + 1*TMath::Pi();
1635   } // this assumes we are doing full jets currently
1636
1637   if( (dphi>0) && (dphi<1*TMath::Pi()/2) ){
1638     // Do nothing! we are in quadrant 1
1639   }else if( (dphi>1*TMath::Pi()/2) && (dphi<1*TMath::Pi()) ){
1640     dphi = 1*TMath::Pi() - dphi;
1641   }else if( (dphi<0) && (dphi>-1*TMath::Pi()/2) ){
1642     dphi = fabs(dphi);
1643   }else if( (dphi<-1*TMath::Pi()/2) && (dphi>-1*TMath::Pi()) ){
1644     dphi = dphi + 1*TMath::Pi();
1645   }
1646
1647   // test
1648   if( dphi < 0 || dphi > TMath::Pi()/2 )
1649     AliWarning(Form("%s: dPHI not in range [0, 0.5*Pi]!", GetName()));
1650
1651   return dphi;   // dphi in [0, Pi/2]
1652 }
1653
1654 //Int_t ieta=GetEtaBin(deta);
1655 //________________________________________________________________________
1656 Int_t AliAnalysisTaskEmcalJetHadEPpid::GetEtaBin(Double_t eta) const
1657 {
1658   // Get eta bin for histos.
1659   Int_t etabin = -1;
1660   if (TMath::Abs(eta)<=0.4)                                             etabin = 0;
1661   else if (TMath::Abs(eta)>0.4 && TMath::Abs(eta)<0.8)  etabin = 1;
1662   else if (TMath::Abs(eta)>=0.8)                                    etabin = 2;
1663
1664   return etabin;
1665 } // end of get-eta-bin
1666
1667 //________________________________________________________________________
1668 Int_t AliAnalysisTaskEmcalJetHadEPpid::GetpTjetBin(Double_t pt) const
1669 {
1670   // Get jet pt  bin for histos.
1671   Int_t ptbin = -1;
1672   if (pt>=15 && pt<20)          ptbin = 0;
1673   else if (pt>=20 && pt<25)     ptbin = 1;
1674   else if (pt>=25 && pt<40)     ptbin = 2;
1675   else if (pt>=40 && pt<60)     ptbin = 3;
1676   else if (pt>=60)              ptbin = 4;
1677
1678   return ptbin;
1679 } // end of get-jet-pt-bin
1680
1681 //________________________________________________________________________
1682 Int_t AliAnalysisTaskEmcalJetHadEPpid::GetpTtrackBin(Double_t pt) const
1683 {
1684   // May need to update bins for future runs... (testing locally)
1685
1686   // Get track pt bin for histos.
1687   Int_t ptbin = -1;
1688   if (pt < 0.5)                 ptbin = 0;
1689   else if (pt>=0.5 && pt<1.0)   ptbin = 1;
1690   else if (pt>=1.0 && pt<1.5)   ptbin = 2;
1691   else if (pt>=1.5 && pt<2.0)   ptbin = 3;
1692   else if (pt>=2.0 && pt<2.5)   ptbin = 4;
1693   else if (pt>=2.5 && pt<3.0)   ptbin = 5;
1694   else if (pt>=3.0 && pt<4.0)   ptbin = 6;
1695   else if (pt>=4.0 && pt<5.0)   ptbin = 7;
1696   else if (pt>=5.0)             ptbin = 8;
1697
1698   return ptbin;
1699 } // end of get-jet-pt-bin
1700
1701 //___________________________________________________________________________
1702 Int_t AliAnalysisTaskEmcalJetHadEPpid::GetzVertexBin(Double_t zVtx) const
1703 {
1704   // get z-vertex bin for histo.
1705   int zVbin= -1;
1706   if (zVtx>=-10 && zVtx<-8)         zVbin = 0;
1707   else if (zVtx>=-8 && zVtx<-6) zVbin = 1;
1708   else if (zVtx>=-6 && zVtx<-4) zVbin = 2;
1709   else if (zVtx>=-4 && zVtx<-2) zVbin = 3;
1710   else if (zVtx>=-2 && zVtx<0)  zVbin = 4;
1711   else if (zVtx>=0 && zVtx<2)   zVbin = 5;
1712   else if (zVtx>=2 && zVtx<4)   zVbin = 6;
1713   else if (zVtx>=4 && zVtx<6)   zVbin = 7;
1714   else if (zVtx>=6 && zVtx<8)   zVbin = 8;
1715   else if (zVtx>=8 && zVtx<10)  zVbin = 9;
1716   else zVbin = 10;
1717
1718   return zVbin;
1719 } // end of get z-vertex bin
1720
1721 //______________________________________________________________________
1722 THnSparse* AliAnalysisTaskEmcalJetHadEPpid::NewTHnSparseF(const char* name, UInt_t entries)
1723 {
1724    // generate new THnSparseF, axes are defined in GetDimParams()
1725    Int_t count = 0;
1726    UInt_t tmp = entries;
1727    while(tmp!=0){
1728       count++;
1729       tmp = tmp &~ -tmp;  // clear lowest bit
1730    }
1731
1732    TString hnTitle(name);
1733    const Int_t dim = count;
1734    Int_t nbins[dim];
1735    Double_t xmin[dim];
1736    Double_t xmax[dim];
1737
1738    Int_t i=0;
1739    Int_t c=0;
1740    while(c<dim && i<32){
1741       if(entries&(1<<i)){
1742
1743          TString label("");
1744          GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
1745          hnTitle += Form(";%s",label.Data());
1746          c++;
1747       }
1748
1749       i++;
1750    }
1751    hnTitle += ";";
1752
1753    return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
1754 } // end of NewTHnSparseF
1755
1756 void AliAnalysisTaskEmcalJetHadEPpid::GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
1757 {
1758    // stores label and binning of axis for THnSparse
1759    const Double_t pi = TMath::Pi();
1760
1761    switch(iEntry){
1762
1763    case 0:
1764       label = "V0 centrality (%)";
1765       nbins = 10;
1766       xmin = 0.;
1767       xmax = 100.;
1768       break;
1769
1770    case 1:
1771       label = "Jet p_{T}";
1772       nbins = 50; // 216
1773       xmin = 0.;
1774       xmax = 250.; // 216
1775       break;
1776
1777    case 2:
1778       label = "Track p_{T}";
1779       nbins = 80; //300; // 750 pid
1780       xmin = 0.;
1781       xmax = 20.; //75.;
1782       break;
1783
1784     case 3:
1785       label = "Relative Eta";
1786       nbins = 48;
1787       xmin = -1.6;
1788       xmax = 1.6;
1789       break;
1790
1791    case 4:
1792       label = "Relative Phi";
1793       nbins = 72;
1794       xmin = -0.5*pi;
1795       xmax = 1.5*pi;
1796       break;
1797
1798   case 5:
1799       label = "Relative angle of Jet and Reaction Plane";
1800       nbins = 12; // 72
1801       xmin = 0;
1802       xmax = 0.5*pi;
1803       break;
1804
1805   case 6:
1806       label = "z-vertex";
1807       nbins = 10;
1808       xmin = -10;
1809       xmax =  10;
1810       break;
1811
1812   case 7:
1813       label = "track charge";
1814       nbins = 3;
1815       xmin = -1.5;
1816       xmax = 1.5;
1817       break;
1818
1819   case 8:
1820       label = "leading jet";
1821       nbins = 3;
1822       xmin = -0.5;
1823       xmax = 2.5;
1824       break;
1825
1826   case 9: // need to update
1827       label = "leading track";
1828       nbins = 10;
1829       xmin = 0;
1830       xmax = 50;
1831       break;
1832
1833    } // end of switch
1834 } // end of getting dim-params
1835
1836 //_________________________________________________
1837 // From CF event mixing code PhiCorrelations
1838 TObjArray* AliAnalysisTaskEmcalJetHadEPpid::CloneAndReduceTrackList(TObjArray* tracksME)
1839 {
1840   // clones a track list by using AliPicoTrack which uses much less memory (used for event mixing)
1841   TObjArray* tracksClone = new TObjArray;
1842   tracksClone->SetOwner(kTRUE);
1843
1844 // ===============================
1845
1846 //      cout << "RM Hybrid track : " << i << "  " << particle->Pt() << endl;
1847
1848   //cout << "nEntries " << tracks->GetEntriesFast() <<endl;
1849   for (Int_t i=0; i<tracksME->GetEntriesFast(); i++) {         // AOD/general case
1850     AliVParticle* particle = (AliVParticle*) tracksME->At(i);  // AOD/general case
1851     if(TMath::Abs(particle->Eta())>fTrkEta) continue;
1852     if(particle->Pt()<0.15)continue;
1853
1854 /*
1855 // DON'T USE
1856     Double_t trackpt=particle->Pt();   // track pT
1857
1858     Int_t trklabel=-1;
1859     trklabel=particle->GetLabel();
1860     //cout << "TRACK_LABEL: " << particle->GetLabel()<<endl;
1861
1862     Int_t hadbin=-1;
1863     if(trackpt<0.5) hadbin=0;
1864     else if(trackpt<1) hadbin=1;
1865     else if(trackpt<2) hadbin=2;
1866     else if(trackpt<3) hadbin=3;
1867     else if(trackpt<5) hadbin=4;
1868     else if(trackpt<8) hadbin=5;
1869     else if(trackpt<20) hadbin=6;
1870 // end of DON'T USE
1871
1872 //feb10 comment out
1873     if(hadbin>-1 && trklabel>-1 && trklabel <3) fHistTrackEtaPhi[trklabel][hadbin]->Fill(particle->Eta(),particle->Phi());
1874     if(hadbin>-1) fHistTrackEtaPhi[3][hadbin]->Fill(particle->Eta(),particle->Phi());
1875
1876     if(hadbin>-1) fHistTrackEtaPhi[hadbin]->Fill(particle->Eta(),particle->Phi());  // TEST
1877 */
1878
1879     tracksClone->Add(new AliPicoTrack(particle->Pt(), particle->Eta(), particle->Phi(), particle->Charge(), 0, 0, 0, 0));
1880   } // end of looping through tracks
1881
1882   return tracksClone;
1883 }
1884
1885 //____________________________________________________________________________________________
1886 THnSparse* AliAnalysisTaskEmcalJetHadEPpid::NewTHnSparseFPID(const char* name, UInt_t entries)
1887 {
1888    // generate new THnSparseF PID, axes are defined in GetDimParams()
1889    Int_t count = 0;
1890    UInt_t tmp = entries;
1891    while(tmp!=0){
1892       count++;
1893       tmp = tmp &~ -tmp;  // clear lowest bit
1894    }
1895
1896    TString hnTitle(name);
1897    const Int_t dim = count;
1898    Int_t nbins[dim];
1899    Double_t xmin[dim];
1900    Double_t xmax[dim];
1901
1902    Int_t i=0;
1903    Int_t c=0;
1904    while(c<dim && i<32){
1905       if(entries&(1<<i)){
1906
1907          TString label("");
1908          GetDimParamsPID(i, label, nbins[c], xmin[c], xmax[c]);
1909          hnTitle += Form(";%s",label.Data());
1910          c++;
1911       }
1912
1913       i++;
1914    }
1915    hnTitle += ";";
1916
1917    return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
1918 } // end of NewTHnSparseF PID
1919
1920 //________________________________________________________________________________
1921 void AliAnalysisTaskEmcalJetHadEPpid::GetDimParamsPID(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
1922 {
1923    // stores label and binning of axis for THnSparse
1924    const Double_t pi = TMath::Pi();
1925
1926    switch(iEntry){
1927
1928    case 0:
1929       label = "V0 centrality (%)";
1930       nbins = 10;
1931       xmin = 0.;
1932       xmax = 100.;
1933       break;
1934
1935    case 1:
1936       label = "Track p_{T}";
1937       nbins = 80; //300; // 750
1938       xmin = 0.;
1939       xmax = 20.; //75.;
1940       break;
1941
1942    case 2:
1943       label = "Charge of Track";
1944       nbins = 3;
1945       xmin = -1.5;
1946       xmax = 1.5;
1947       break;
1948
1949    case 3:
1950       label = "Relative Eta of Track and Jet";
1951       nbins = 48;
1952       xmin = -1.6;
1953       xmax = 1.6;
1954       break;
1955
1956    case 4:
1957       label = "Relative Phi of Track and Jet";
1958       nbins = 72;
1959       xmin = -0.5*pi;
1960       xmax = 1.5*pi;
1961       break;
1962
1963    case 5:
1964       label = "leading jet";
1965       nbins = 3;
1966       xmin = -.5;
1967       xmax = 2.5;
1968       break;
1969
1970    case 6:
1971       label = "Z-vertex";
1972       nbins = 10;
1973       xmin = -10.;
1974       xmax = 10.;
1975       break;
1976
1977    case 7:
1978       label = "Relative angle: Jet and Reaction Plane";
1979       nbins = 12; // 48
1980       xmin = 0.;
1981       xmax = 0.5*pi;
1982       break;
1983
1984    case 8:
1985       label = "Jet p_{T}";
1986       nbins = 50; // 216
1987       xmin = 0.;
1988       xmax = 250.; // 216
1989       break;
1990
1991    case 9:
1992       label = "N-Sigma of pions in TPC";
1993       nbins = 200;
1994       xmin = -10.0;
1995       xmax = 10.0;
1996       break;
1997
1998    case 10:
1999       label = "N-Sigma of pions in TOF";
2000       nbins = 200;
2001       xmin = -10.;
2002       xmax = 10.;
2003       break;
2004
2005    case 11:
2006       label = "TPC PID determination";
2007       nbins = 5;
2008       xmin = 0.;
2009       xmax = 5.;
2010       break;
2011
2012    case 12:
2013       label = "ITS PID determination";
2014       nbins = 5;
2015       xmin = 0.;
2016       xmax = 5.;
2017       break;
2018
2019    case 13:
2020       label = "TOF PID determination";
2021       nbins = 5;
2022       xmin = 0.;
2023       xmax = 5.;
2024       break;
2025
2026    case 14:
2027       label = "N-Sigma of protons in TPC";
2028       nbins = 200;
2029       xmin = -10.;
2030       xmax = 10.;
2031       break;
2032
2033    case 15:
2034       label = "N-Sigma of kaons in TPC";
2035       nbins = 200;
2036       xmin = -10.;
2037       xmax = 10.;
2038       break;
2039
2040    case 16:
2041       label = "N-Sigma of pions in ITS";
2042       nbins = 200;
2043       xmin = -10.0;
2044       xmax = 10.0;
2045       break;
2046
2047    case 17:
2048       label = "N-Sigma of protons in ITS";
2049       nbins = 200;
2050       xmin = -10.;
2051       xmax = 10.;
2052       break;
2053
2054    case 18:
2055       label = "N-Sigma of kaons in ITS";
2056       nbins = 200;
2057       xmin = -10.;
2058       xmax = 10.;
2059       break;
2060
2061    case 19:
2062       label = "N-Sigma of protons in TOF";
2063       nbins = 200;
2064       xmin = -10.;
2065       xmax = 10.;
2066       break;
2067
2068    case 20:
2069       label = "N-Sigma of kaons in TOF";
2070       nbins = 200;
2071       xmin = -10.;
2072       xmax = 10.;
2073       break;
2074
2075    } // end of switch
2076 } // end of get dimension parameters PID
2077
2078 void AliAnalysisTaskEmcalJetHadEPpid::Terminate(Option_t *) {
2079   cout<<"#########################"<<endl;
2080   cout<<"#### DONE RUNNING!!! ####"<<endl;
2081   cout<<"#########################"<<endl;
2082 } // end of terminate
2083
2084 //________________________________________________________________________
2085 Int_t AliAnalysisTaskEmcalJetHadEPpid::AcceptMyJet(AliEmcalJet *jet) {
2086   //applies all jet cuts except pt
2087   if ((jet->Phi()<fPhimin)||(jet->Phi()>fPhimax)) return 0;
2088   if ((jet->Eta()<fEtamin)||(jet->Eta()>fEtamax)) return 0;
2089   if (jet->Area()<fAreacut) return 0;
2090   // prevents 0 area jets from sneaking by when area cut == 0
2091   if (jet->Area()==0) return 0;
2092   //exclude jets with extremely high pt tracks which are likely misreconstructed
2093   if(jet->MaxTrackPt()>100) return 0;
2094
2095   //passed all above cuts
2096   return 1;
2097 }
2098
2099 //void AliAnalysisTaskEmcalJetHadEPpid::FillAnalysisSummaryHistogram() const
2100 void AliAnalysisTaskEmcalJetHadEPpid::SetfHistPIDcounterLabels(TH1* h) const
2101 {
2102     // fill the analysis summary histrogram, saves all relevant analysis settigns
2103     h->GetXaxis()->SetBinLabel(1, "TPC: Unidentified");
2104     h->GetXaxis()->SetBinLabel(2, "TPC: Pion");
2105     h->GetXaxis()->SetBinLabel(3, "TPC: Kaon");
2106     h->GetXaxis()->SetBinLabel(4, "TPC: Proton");
2107     h->GetXaxis()->SetBinLabel(5, "ITS: Unidentified");
2108     h->GetXaxis()->SetBinLabel(6, "ITS: Pion");
2109     h->GetXaxis()->SetBinLabel(7, "ITS: Kaon");
2110     h->GetXaxis()->SetBinLabel(8, "ITS: Proton");
2111     h->GetXaxis()->SetBinLabel(9, "TOF: Unidentified");
2112     h->GetXaxis()->SetBinLabel(10, "TOF: Pion");
2113     h->GetXaxis()->SetBinLabel(11, "TOF: Kaon");
2114     h->GetXaxis()->SetBinLabel(12, "TOF: Proton");
2115     h->GetXaxis()->SetBinLabel(14, "Unidentified tracks");
2116
2117     // set x-axis labels vertically
2118     h->LabelsOption("v");
2119 }
2120
2121 //void AliAnalysisTaskEmcalJetHadEPpid::FillAnalysisSummaryHistogram() const
2122 void AliAnalysisTaskEmcalJetHadEPpid::SetfHistQAcounterLabels(TH1* h) const
2123 {
2124     // label bins of the analysis event summary
2125     h->GetXaxis()->SetBinLabel(1, "All events started");
2126     h->GetXaxis()->SetBinLabel(2, "object check");
2127     h->GetXaxis()->SetBinLabel(3, "aod/esd check");
2128     h->GetXaxis()->SetBinLabel(4, "centrality check");
2129     h->GetXaxis()->SetBinLabel(5, "zvertex check");
2130     h->GetXaxis()->SetBinLabel(6, "list check");
2131     h->GetXaxis()->SetBinLabel(7, "track/jet pointer check");
2132     h->GetXaxis()->SetBinLabel(8, "tracks & jets < than 1 check");
2133     h->GetXaxis()->SetBinLabel(9, "after track/jet loop to get highest pt");
2134     h->GetXaxis()->SetBinLabel(10, "accepted jets");
2135     h->GetXaxis()->SetBinLabel(11, "jets meeting pt threshold");
2136     h->GetXaxis()->SetBinLabel(12, "accepted tracks in events w/ trigger jet");
2137     h->GetXaxis()->SetBinLabel(13, "after AliVEvent & fPIDResponse");
2138     h->GetXaxis()->SetBinLabel(14, "events before event mixing");
2139     h->GetXaxis()->SetBinLabel(15, "mixed events w/ pool");
2140     h->GetXaxis()->SetBinLabel(16, "event mixing: jets");
2141     h->GetXaxis()->SetBinLabel(17, "event mixing: nMix");
2142     h->GetXaxis()->SetBinLabel(18, "event mixing: nbackground tracks");
2143     h->GetXaxis()->SetBinLabel(19, "event mixing: THE END");
2144
2145     // set x-axis labels vertically
2146     h->LabelsOption("v");
2147 }
2148
2149 void AliAnalysisTaskEmcalJetHadEPpid::SetfHistEvtSelQALabels(TH1* h) const
2150 {
2151     // label bins of the analysis trigger selection summary
2152     h->GetXaxis()->SetBinLabel(1, "no trigger");
2153     h->GetXaxis()->SetBinLabel(2, "kAny");
2154     h->GetXaxis()->SetBinLabel(3, "kAnyINT");
2155     h->GetXaxis()->SetBinLabel(4, "kMB");
2156     h->GetXaxis()->SetBinLabel(5, "kINT7");
2157     h->GetXaxis()->SetBinLabel(6, "kEMC1");
2158     h->GetXaxis()->SetBinLabel(7, "kEMC7");
2159     h->GetXaxis()->SetBinLabel(8, "kEMC8");
2160     h->GetXaxis()->SetBinLabel(9, "kEMCEJE");
2161     h->GetXaxis()->SetBinLabel(10, "kEMCEGA");
2162     h->GetXaxis()->SetBinLabel(11, "kCentral");
2163     h->GetXaxis()->SetBinLabel(12, "kSemiCentral");
2164     h->GetXaxis()->SetBinLabel(13, "kINT8");
2165     h->GetXaxis()->SetBinLabel(14, "kEMCEJE or kMB");
2166     h->GetXaxis()->SetBinLabel(15, "kEMCEGA or kMB");
2167     h->GetXaxis()->SetBinLabel(16, "kAnyINT or kMB");
2168     h->GetXaxis()->SetBinLabel(17, "kEMCEJE & kMB");
2169     h->GetXaxis()->SetBinLabel(18, "kEMCEGA & kMB");
2170     h->GetXaxis()->SetBinLabel(19, "kAnyINT & kMB");
2171
2172
2173     // set x-axis labels vertically
2174     h->LabelsOption("v");
2175     //h->LabelsDeflate("X");
2176 }
2177
2178 //______________________________________________________________________
2179 THnSparse* AliAnalysisTaskEmcalJetHadEPpid::NewTHnSparseFCorr(const char* name, UInt_t entries) {
2180   // generate new THnSparseD, axes are defined in GetDimParamsD()
2181   Int_t count = 0;
2182   UInt_t tmp = entries;
2183   while(tmp!=0){
2184     count++;
2185     tmp = tmp &~ -tmp;  // clear lowest bit
2186   }
2187
2188   TString hnTitle(name);
2189   const Int_t dim = count;
2190   Int_t nbins[dim];
2191   Double_t xmin[dim];
2192   Double_t xmax[dim];
2193
2194   Int_t i=0;
2195   Int_t c=0;
2196   while(c<dim && i<32){
2197     if(entries&(1<<i)){
2198
2199       TString label("");
2200       GetDimParamsCorr(i, label, nbins[c], xmin[c], xmax[c]);
2201       hnTitle += Form(";%s",label.Data());
2202       c++;
2203     }
2204
2205     i++;
2206   }
2207   hnTitle += ";";
2208
2209   return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
2210 } // end of NewTHnSparseF
2211
2212 //______________________________________________________________________________________________
2213 void AliAnalysisTaskEmcalJetHadEPpid::GetDimParamsCorr(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
2214 {
2215    //stores label and binning of axis for THnSparse
2216    const Double_t pi = TMath::Pi();
2217
2218    switch(iEntry){
2219
2220     case 0:
2221       label = "V0 centrality (%)";
2222       nbins = 10;
2223       xmin = 0.;
2224       xmax = 100.;
2225       break;
2226
2227    case 1:
2228       label = "Jet p_{T}";
2229       nbins = 50; // 216
2230       xmin = 0.;
2231       xmax = 250.; // 216
2232       break;
2233
2234    case 2:
2235       label = "Relative angle: Jet and Reaction Plane";
2236       nbins = 12; // 48
2237       xmin = 0.;
2238       xmax = 0.5*pi;
2239       break;
2240
2241    case 3:
2242           label = "Z-vertex";
2243       nbins = 10;
2244       xmin = -10.;
2245       xmax = 10.;
2246           break;
2247
2248    case 4:
2249           label = "Jet p_{T} corrected with Local Rho";
2250       nbins = 50; // 250
2251       xmin = -50.;
2252       xmax = 200.;
2253           break;
2254
2255    case 5:
2256           label = "Jet p_{T} corrected with Global Rho";
2257       nbins = 50; // 250
2258       xmin = -50.;
2259       xmax = 200.;
2260           break;
2261
2262    }// end of switch
2263 } // end of Correction (ME) sparse
2264
2265 //________________________________________________________________________
2266 //Int_t AliAnalysisTaskEmcalJetHadEPpid::AcceptFlavourJet(AliEmcalJet* fljet, Int_t NUM, Int_t NUM2, Int_t NUM3) {
2267 Int_t AliAnalysisTaskEmcalJetHadEPpid::AcceptFlavourJet(AliEmcalJet* fljet, Int_t NUM) {
2268   // Get jet if accepted for given flavour tag
2269   // If jet not accepted return 0
2270   if(!fljet) {
2271     AliError(Form("%s:Jet not found",GetName()));
2272     return 0;
2273   }
2274
2275   Int_t flavNUM = -99;//, flavNUM2 = -99, flavNUM3 = -99; FIXME commented out to avoid compiler warning
2276   flavNUM = NUM;
2277   //flavNUM2 = NUM2;
2278   //flavNUM3 = NUM3;
2279
2280 /*
2281   // from the AliEmcalJet class, the tagging enumerator
2282   enum EFlavourTag{
2283        kDStar = 1<<0, kD0 = 1<<1,
2284            kSig1 = 1<<2, kSig2 = 1<<3,
2285            kBckgrd1 = 1<<4, kBckgrd2 = 1<<5, kBckgrd3 = 1<<6
2286   };
2287   // bit 0 = no tag, bit 1 = Dstar, bit 2 = D0 and so forth...
2288 */
2289
2290   // get flavour of jet (if any)
2291   Int_t flav = -999;
2292   flav = fljet->GetFlavour();
2293
2294   // cases (for now..)
2295   // 3 = electron rich, 5 = hadron (bkgrd) rich
2296   // if flav < 1, its not tagged, so return kFALSE (0)
2297   if(flav < 1) return 0;
2298
2299   // if flav is not equal to what we want then return kFALSE (0)
2300   //if(flav != flavNUM && flav != flavNUM2 && flav != flavNUM3) return 0;
2301   if(flav != flavNUM) return 0;
2302
2303   // we have the flavour we want, so return kTRUE (1)
2304   //if(flav == flavNUM || flav == flavNUM2 || flav == flavNUM3) return 1;
2305   if(flav == flavNUM) return 1;
2306
2307   // we by default have a flavour tagged jet
2308   return 1;
2309 }