PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetHadEPpid.cxx

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