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