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