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