AliPhysics  8bb951a (8bb951a)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskEmcalJetHadEPpid.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 //-------------------------------------------------------------------------
17 // 1) Analysis Task to perform Jet-Hadron Correlations
18 // 2) Event plane dependence task.
19 // 3) performs event plane resolution calculation
20 // 4) does PID of the associated pi/k/p hadrons
21 //
22 // Author: Joel Mazer (joel.mazer@cern.ch)
23 //-------------------------------------------------------------------------
24 
25 // task head include
27 
28 // general ROOT includes
29 #include <TCanvas.h>
30 #include <TMath.h>
31 #include <TChain.h>
32 #include <TClonesArray.h>
33 #include <TH1F.h>
34 #include <TH2F.h>
35 #include <TH3F.h>
36 #include <TH1.h>
37 #include <TH2.h>
38 #include <TH3.h>
39 #include <TProfile.h>
40 #include <THnSparse.h>
41 #include <TList.h>
42 #include <TLorentzVector.h>
43 #include <TParameter.h>
44 #include <TParticle.h>
45 #include <TTree.h>
46 #include <TVector3.h>
47 #include <TObjArray.h>
48 
49 // AliROOT includes
50 #include "AliAODEvent.h"
51 #include "AliESDEvent.h"
52 #include "AliAnalysisManager.h"
53 #include "AliAnalysisTask.h"
54 #include "AliCentrality.h"
55 #include "AliEmcalJet.h"
56 #include "AliAODJet.h"
57 #include "AliVCluster.h"
58 #include "AliVTrack.h"
59 #include <AliVEvent.h>
60 #include <AliVParticle.h>
61 #include "AliRhoParameter.h"
62 #include "AliEmcalParticle.h"
63 
64 // Localized Rho includes
65 #include "AliLocalRhoParameter.h"
67 
68 // event handler (and pico's) includes
69 #include <AliInputEventHandler.h>
70 #include <AliVEventHandler.h>
71 #include "AliESDInputHandler.h"
72 #include "AliESDCaloCluster.h"
73 #include "AliPicoTrack.h"
74 #include "AliEventPoolManager.h"
75 #include "AliESDtrackCuts.h"
76 
77 // PID includes
78 #include "AliPIDResponse.h"
79 #include "AliTPCPIDResponse.h"
80 #include "AliESDpid.h"
81 
82 // magnetic field includes
83 #include "TGeoGlobalMagField.h"
84 #include "AliMagF.h"
85 
86 // container includes
87 #include "AliJetContainer.h"
88 #include "AliParticleContainer.h"
89 #include "AliClusterContainer.h"
90 
91 using std::cout;
92 using std::endl;
93 
95 
96 //________________________________________________________________________
98  AliAnalysisTaskEmcalJet("correlations",kFALSE),
99  fPhimin(-10), fPhimax(10),
100  fEtamin(-0.9), fEtamax(0.9),
101  fAreacut(0.0), fTrkBias(5), fClusBias(5), fTrkEta(0.9),
102  fJetPtcut(15.0), fJetRad(0.4), fConstituentCut(3.0),
103  fesdTrackCuts(0),
104  fDetectorType(kVZEROComb),
105  fSoftTrackMinPt_ep(0.15), fSoftTrackMaxPt_ep(5.),
106  fNAcceptedTracks(0), fLeadingJet(0), fExcludeLeadingJetsFromFit(1.),
107  fCentralityClasses(0), fInCentralitySelection(-1),
108  fDoEventMixing(0), fMixingTracks(50000), fNMIXtracks(5000), fNMIXevents(5),
109  fCentBinSize(1), fReduceStatsCent(0),
110  fTriggerEventType(AliVEvent::kAny), fMixingEventType(AliVEvent::kAny),
111  fDoEffCorr(0), fcorrJetPt(0),
112  doPlotGlobalRho(0), doVariableBinning(0), dovarbinTHnSparse(0),
113  makeQAhistos(0), makeBIAShistos(0), makeextraCORRhistos(0), makeoldJEThadhistos(0),
114  allpidAXIS(0), fcutType("EMCAL"), doPID(0), doPIDtrackBIAS(0), doaltPIDbinning(0),
115  doEventPlaneRes(0),
116  doComments(0),
117  doFlavourJetAnalysis(0), fJetFlavTag(3),
118  fBeam(kNA),
119  fLocalRhoVal(0),
120  fTracksName(""), fTracksNameME(""), fJetsName(""), fCaloClustersName(""),
121  event(0),
122  isPItpc(0), isKtpc(0), isPtpc(0), // pid TPC
123  isPIits(0), isKits(0), isPits(0), // pid ITS
124  isPItof(0), isKtof(0), isPtof(0), // pid TOF
125  fPoolMgr(0x0),
126  fPIDResponse(0x0), fTPCResponse(),
127  fESD(0), fAOD(0), fVevent(0),
128  fHistEventQA(0), fHistEventSelectionQA(0), fHistEventSelectionQAafterCuts(0),
129  fHistCentZvertGA(0), fHistCentZvertJE(0), fHistCentZvertMB(0), fHistCentZvertAny(0),
130  fHistTPCdEdX(0), fHistITSsignal(0), //fHistTOFsignal(0),
131  fHistRhovsCent(0), fHistNjetvsCent(0),
132  fHistMult(0),
133  fHistJetPhi(0), fHistTrackPhi(0),
134  fHistLocalRhoJetpt(0),
135  fHistJetHaddPhiIN(0), fHistJetHaddPhiOUT(0), fHistJetHaddPhiMID(0),
136  fHistJetHaddPhiBias(0), fHistJetHaddPhiINBias(0), fHistJetHaddPhiOUTBias(0), fHistJetHaddPhiMIDBias(0),
137  fHistMEdPHI(0),
138  fHistTrackPtallcent(0),
139  fHistJetEtaPhi(0), fHistJetHEtaPhi(0),
140  fHistSEphieta(0), fHistMEphieta(0),
141  fHistJetHaddPHI(0),
142  fHistClusEtaPhiEnergy(0), fHistClusEnergy(0), fHistClusofJetEnergy(0),
143  fHistJetNClusterConstit(0), fHistJetNTrackConstit(0), fHistJetNConstit(0),
144  fHistPID(0),
145  fhnPID(0x0), fhnMixedEvents(0x0), fhnJH(0x0), fhnCorr(0x0),
146  fChi2A(0x0), fChi2C(0x0), fChi3A(0x0), fChi3C(0x0), fUseChiWeightForVZERO(kTRUE), // test
147  fTracksFromContainer(0),
148  fJetsCont(0), fTracksCont(0), fCaloClustersCont(0),
149  fContainerAllJets(0), fContainerPIDJets(1)//, fTrgJet(0)
150 {
151  // Default Constructor
152  for(Int_t ilab=0; ilab<4; ilab++){
153  for(Int_t ipta=0; ipta<7; ipta++){
154  //fHistTrackEtaPhi[ilab][ipta]=0; // keep out for now
155  } // end of pt-associated loop
156  } // end of lab loop
157 
158  for(Int_t cen = 0; cen<10; cen++){
159  fProfV2Resolution[cen]=0;
160  fProfV3Resolution[cen]=0;
161  fProfV4Resolution[cen]=0;
162  fProfV5Resolution[cen]=0;
163  }
164 
165  for(Int_t itrackpt=0; itrackpt<9; itrackpt++){
166  fHistJetHadbindPhi[itrackpt]=0;
167  fHistJetHadbindPhiIN[itrackpt]=0;
168  fHistJetHadbindPhiMID[itrackpt]=0;
169  fHistJetHadbindPhiOUT[itrackpt]=0;
170  } // end of trackpt bin loop
171 
172  for(Int_t icent = 0; icent<6; ++icent){
173  for(Int_t iptjet = 0; iptjet<5; ++iptjet){
174  for(Int_t ieta = 0; ieta<3; ++ieta){
175  fHistJetH[icent][iptjet][ieta]=0;
176  fHistJetHBias[icent][iptjet][ieta]=0;
177  fHistJetHTT[icent][iptjet][ieta]=0;
178  } // end of eta loop
179  } // end of pt-jet loop
180  } // end of centrality loop
181 
182  // centrality dependent histo's
183  for (Int_t i = 0;i<6;++i){
184  fHistJetPt[i] = 0;
185  fHistJetPtBias[i] = 0;
186  fHistJetPtTT[i] = 0;
187  fHistAreavsRawPt[i] = 0;
188  fHistJetPtvsTrackPt[i] = 0;
189  fHistTrackPt[i] = 0;
190  fHistEP0[i] = 0;
191  fHistEP0A[i] = 0;
192  fHistEP0C[i] = 0;
193  fHistEPAvsC[i] = 0;
194  fHistJetPtcorrGlRho[i] = 0;
195  fHistJetPtvsdEP[i] = 0;
196  fHistJetPtvsdEPBias[i] = 0;
197  fHistRhovsdEP[i] = 0;
198  fHistJetEtaPhiPt[i] = 0;
199  fHistJetEtaPhiPtBias[i] = 0;
200  fHistJetPtArea[i] = 0;
201  fHistJetPtAreaBias[i] = 0;
202  fHistJetPtNcon[i] = 0;
203  fHistJetPtNconBias[i] = 0;
204  fHistJetPtNconCh[i] = 0;
205  fHistJetPtNconBiasCh[i] = 0;
206  fHistJetPtNconEm[i] = 0;
207  fHistJetPtNconBiasEm[i] = 0;
208  fHistJetHaddPhiINcent[i] = 0;
209  fHistJetHaddPhiOUTcent[i] = 0;
210  fHistJetHaddPhiMIDcent[i] = 0;
211  }
212 
213  SetMakeGeneralHistograms(kTRUE);
214 
215 }
216 
217 //________________________________________________________________________
219  AliAnalysisTaskEmcalJet(name,kTRUE),
220  fPhimin(-10), fPhimax(10),
221  fEtamin(-0.9), fEtamax(0.9),
222  fAreacut(0.0), fTrkBias(5), fClusBias(5), fTrkEta(0.9),
223  fJetPtcut(15.0), fJetRad(0.4), fConstituentCut(3.0),
224  fesdTrackCuts(0),
225  fDetectorType(kVZEROComb),
226  fSoftTrackMinPt_ep(0.15), fSoftTrackMaxPt_ep(5.),
227  fNAcceptedTracks(0), fLeadingJet(0), fExcludeLeadingJetsFromFit(1.),
228  fCentralityClasses(0), fInCentralitySelection(-1),
229  fDoEventMixing(0), fMixingTracks(50000), fNMIXtracks(5000), fNMIXevents(5),
230  fCentBinSize(1), fReduceStatsCent(0),
231  fTriggerEventType(AliVEvent::kAny), fMixingEventType(AliVEvent::kAny),
232  fDoEffCorr(0), fcorrJetPt(0),
233  doPlotGlobalRho(0), doVariableBinning(0), dovarbinTHnSparse(0),
234  makeQAhistos(0), makeBIAShistos(0), makeextraCORRhistos(0), makeoldJEThadhistos(0),
235  allpidAXIS(0), fcutType("EMCAL"), doPID(0), doPIDtrackBIAS(0), doaltPIDbinning(0),
236  doEventPlaneRes(0),
237  doComments(0),
238  doFlavourJetAnalysis(0), fJetFlavTag(3),
239  fBeam(kNA),
240  fLocalRhoVal(0),
241  fTracksName(""), fTracksNameME(""), fJetsName(""), fCaloClustersName(""),
242  event(0),
243  isPItpc(0), isKtpc(0), isPtpc(0), // pid TPC
244  isPIits(0), isKits(0), isPits(0), // pid ITS
245  isPItof(0), isKtof(0), isPtof(0), // pid TOF
246  fPoolMgr(0x0),
247  fPIDResponse(0x0), fTPCResponse(),
248  fESD(0), fAOD(0), fVevent(0),
249  fHistEventQA(0), fHistEventSelectionQA(0), fHistEventSelectionQAafterCuts(0),
250  fHistCentZvertGA(0), fHistCentZvertJE(0), fHistCentZvertMB(0), fHistCentZvertAny(0),
251  fHistTPCdEdX(0), fHistITSsignal(0), //fHistTOFsignal(0),
252  fHistRhovsCent(0), fHistNjetvsCent(0),
253  fHistMult(0),
254  fHistJetPhi(0), fHistTrackPhi(0),
255  fHistLocalRhoJetpt(0),
256  fHistJetHaddPhiIN(0), fHistJetHaddPhiOUT(0), fHistJetHaddPhiMID(0),
257  fHistJetHaddPhiBias(0), fHistJetHaddPhiINBias(0), fHistJetHaddPhiOUTBias(0), fHistJetHaddPhiMIDBias(0),
258  fHistMEdPHI(0),
259  fHistTrackPtallcent(0),
260  fHistJetEtaPhi(0), fHistJetHEtaPhi(0),
261  fHistSEphieta(0), fHistMEphieta(0),
262  fHistJetHaddPHI(0),
263  fHistClusEtaPhiEnergy(0), fHistClusEnergy(0), fHistClusofJetEnergy(0),
264  fHistJetNClusterConstit(0), fHistJetNTrackConstit(0), fHistJetNConstit(0),
265  fHistPID(0),
266  fhnPID(0x0), fhnMixedEvents(0x0), fhnJH(0x0), fhnCorr(0x0),
267  fChi2A(0x0), fChi2C(0x0), fChi3A(0x0), fChi3C(0x0), fUseChiWeightForVZERO(kTRUE), // test
268  fTracksFromContainer(0),
269  fJetsCont(0), fTracksCont(0), fCaloClustersCont(0),
270  fContainerAllJets(0), fContainerPIDJets(1)//, fTrgJet(0)
271 {
272  // Default Constructor
273  for(Int_t ilab=0; ilab<4; ilab++){
274  for(Int_t ipta=0; ipta<7; ipta++){
275  //fHistTrackEtaPhi[ilab][ipta]=0; //keep out for now
276  } // end of pt-associated loop
277  } // end of lab loop
278 
279  for(Int_t cen = 0; cen<10; cen++){
280  fProfV2Resolution[cen]=0;
281  fProfV3Resolution[cen]=0;
282  fProfV4Resolution[cen]=0;
283  fProfV5Resolution[cen]=0;
284  }
285 
286  for(Int_t itrackpt=0; itrackpt<9; itrackpt++){
287  fHistJetHadbindPhi[itrackpt]=0;
288  fHistJetHadbindPhiIN[itrackpt]=0;
289  fHistJetHadbindPhiMID[itrackpt]=0;
290  fHistJetHadbindPhiOUT[itrackpt]=0;
291  } // end of trackpt bin loop
292 
293  for(Int_t icent = 0; icent<6; ++icent){
294  for(Int_t iptjet = 0; iptjet<5; ++iptjet){
295  for(Int_t ieta = 0; ieta<3; ++ieta){
296  fHistJetH[icent][iptjet][ieta]=0;
297  fHistJetHBias[icent][iptjet][ieta]=0;
298  fHistJetHTT[icent][iptjet][ieta]=0;
299  } // end of eta loop
300  } // end of pt-jet loop
301  } // end of centrality loop
302 
303  // centrality dependent histo's
304  for (Int_t i = 0;i<6;++i){
305  fHistJetPt[i] = 0;
306  fHistJetPtBias[i] = 0;
307  fHistJetPtTT[i] = 0;
308  fHistAreavsRawPt[i] = 0;
309  fHistJetPtvsTrackPt[i] = 0;
310  fHistTrackPt[i] = 0;
311  fHistEP0[i] = 0;
312  fHistEP0A[i] = 0;
313  fHistEP0C[i] = 0;
314  fHistEPAvsC[i] = 0;
315  fHistJetPtcorrGlRho[i] = 0;
316  fHistJetPtvsdEP[i] = 0;
317  fHistJetPtvsdEPBias[i] = 0;
318  fHistRhovsdEP[i] = 0;
319  fHistJetEtaPhiPt[i] = 0;
320  fHistJetEtaPhiPtBias[i] = 0;
321  fHistJetPtArea[i] = 0;
322  fHistJetPtAreaBias[i] = 0;
323  fHistJetPtNcon[i] = 0;
324  fHistJetPtNconBias[i] = 0;
325  fHistJetPtNconCh[i] = 0;
326  fHistJetPtNconBiasCh[i] = 0;
327  fHistJetPtNconEm[i] = 0;
328  fHistJetPtNconBiasEm[i] = 0;
329  fHistJetHaddPhiINcent[i] = 0;
330  fHistJetHaddPhiOUTcent[i] = 0;
331  fHistJetHaddPhiMIDcent[i] = 0;
332  }
333 
335 
336 }
337 
338 //_______________________________________________________________________
340 {
341  // destructor
342  if(fOutput) {delete fOutput; fOutput = 0;}
343  if(fCentralityClasses) {delete fCentralityClasses; fCentralityClasses = 0x0;} // test
344  if(fChi2A) {delete fChi2A; fChi2A = 0x0;}
345  if(fChi2C) {delete fChi2C; fChi2C = 0x0;}
346  if(fChi3A) {delete fChi3A; fChi3A = 0x0;}
347  if(fChi3C) {delete fChi3C; fChi3C = 0x0;}
348 }
349 
350 //________________________________________________________________________
352 { // This is called ONCE!
353  if (!fCreateHisto) return;
355  OpenFile(1); // do I need the 1?
356 
357  // char array for naming histograms
358  int nchar = 200;
359  char name[nchar];
360 
361  // strings for titles
362  TString name1;
363  TString title1;
364 
365  // create histograms that arn't array
366  fHistNjetvsCent = new TH2F("NjetvsCent", "NjetvsCent", 100, 0.0, 100.0, 100, 0, 100);
367  fHistJetHaddPHI = new TH1F("fHistJetHaddPHI", "Jet-Hadron #Delta#varphi", 72,-0.5*TMath::Pi(),1.5*TMath::Pi());
368  fHistJetHaddPhiIN = new TH1F("fHistJetHaddPhiIN","Jet-Hadron #Delta#varphi IN PLANE", 72,-0.5*TMath::Pi(), 1.5*TMath::Pi());
369  fHistJetHaddPhiOUT = new TH1F("fHistJetHaddPhiOUT","Jet-Hadron #Delta#varphi OUT PLANE", 72,-0.5*TMath::Pi(), 1.5*TMath::Pi());
370  fHistJetHaddPhiMID = new TH1F("fHistJetHaddPhiMID","Jet-Hadron #Delta#varphi MIDDLE of PLANE", 72,-0.5*TMath::Pi(), 1.5*TMath::Pi());
371  fHistLocalRhoJetpt = new TH1F("fHistLocalRhoJetpt","Local Rho corrected Jet p_{T}", 50, -50, 200);
372 
373  // Centrality and Zvertex distribution for various triggers - Event Mixing QA
374  fHistCentZvertGA = new TH2F("fHistCentZvertGA", "Centrality - Z-vertex distribution for GA trigger", 20, 0, 100, 10, -10, 10);
376  fHistCentZvertJE = new TH2F("fHistCentZvertJE", "Centrality - Z-vertex distribution for JE trigger", 20, 0, 100, 10, -10, 10);
378  fHistCentZvertMB = new TH2F("fHistCentZvertMB", "Centrality - Z-vertex distribution for MB trigger", 20, 0, 100, 10, -10, 10);
380  fHistCentZvertAny = new TH2F("fHistCentZvertAny", "Centrality - Z-vertex distribution for kAny trigger", 20, 0, 100, 10, -10, 10);
382 
383  // Event QA histo
384  fHistEventQA = new TH1F("fHistEventQA", "Event Counter at checkpoints in code", 20, 0.5, 20.5);
386  fOutput->Add(fHistEventQA);
387 
388  // Event Selection QA histo
389  fHistEventSelectionQA = new TH1F("fHistEventSelectionQA", "Trigger Selection Counter", 20, 0.5, 20.5);
390  //SetfHistEvtSelQALabels(fHistEventSelectionQA);
392 
393  // Event Selection QA histo after all cuts to return from task
394  fHistEventSelectionQAafterCuts = new TH1F("fHistEventSelectionQAafterCuts", "Trigger Selection Counter after Cuts", 20, 0.5, 20.5);
395  //SetfHistEvtSelQALabels(fHistEventSelectionQAafterCuts);
397 
398  // add to output lists
399  fOutput->Add(fHistNjetvsCent);
400  fOutput->Add(fHistJetHaddPHI);
405 
406  fHistTPCdEdX = new TH2F("TPCdEdX", "TPCdEdX", 400, 0.0, 20.0, 500, 0, 500);
407  fOutput->Add(fHistTPCdEdX);
408 
409  if(doEventPlaneRes){
410  // Reaction Plane resolution as function of centrality - corrected for 2nd order event plane
411  for (Int_t i = 0; i<10; ++i){
412  fProfV2Resolution[i] = new TProfile(Form("fProfV2Resolution_%i", i), Form("fProfV2Resolution_%i", i), 11, -0.5, 10.5);
413  fProfV2Resolution[i]->GetXaxis()->SetBinLabel(3, "<cos(2(#Psi_{VZEROA} - #Psi_{VZEROC}))>");
414  fProfV2Resolution[i]->GetXaxis()->SetBinLabel(4, "<cos(2(#Psi_{VZEROC} - #Psi_{VZEROA}))>");
415  fProfV2Resolution[i]->GetXaxis()->SetBinLabel(5, "<cos(2(#Psi_{VZEROA} - #Psi_{TPC}))>");
416  fProfV2Resolution[i]->GetXaxis()->SetBinLabel(6, "<cos(2(#Psi_{TPC} - #Psi_{VZEROA}))>");
417  fProfV2Resolution[i]->GetXaxis()->SetBinLabel(7, "<cos(2(#Psi_{VZEROC} - #Psi_{TPC}))>");
418  fProfV2Resolution[i]->GetXaxis()->SetBinLabel(8, "<cos(2(#Psi_{TPC} - #Psi_{VZEROC}))>");
419  fProfV2Resolution[i]->GetXaxis()->SetBinLabel(9, "<cos(2(#Psi_{VZERO} - #Psi_{TPC_A}))>");
420  fProfV2Resolution[i]->GetXaxis()->SetBinLabel(10, "<cos(2(#Psi_{VZERO} - #Psi_{TPC_B}))>");
421  fProfV2Resolution[i]->GetXaxis()->SetBinLabel(11, "<cos(2(#Psi_{TPC_A} - #Psi_{TPC_B}))>");
422  fOutput->Add(fProfV2Resolution[i]);
423  fProfV3Resolution[i] = new TProfile(Form("fProfV3Resolution_%i", i), Form("fProfV3Resolution_%i", i), 11, -0.5, 10.5);
424  fProfV3Resolution[i]->GetXaxis()->SetBinLabel(3, "<cos(3(#Psi_{VZEROA} - #Psi_{VZEROC}))>");
425  fProfV3Resolution[i]->GetXaxis()->SetBinLabel(4, "<cos(3(#Psi_{VZEROC} - #Psi_{VZEROA}))>");
426  fProfV3Resolution[i]->GetXaxis()->SetBinLabel(5, "<cos(3(#Psi_{VZEROA} - #Psi_{TPC}))>");
427  fProfV3Resolution[i]->GetXaxis()->SetBinLabel(6, "<cos(3(#Psi_{TPC} - #Psi_{VZEROA}))>");
428  fProfV3Resolution[i]->GetXaxis()->SetBinLabel(7, "<cos(3(#Psi_{VZEROC} - #Psi_{TPC}))>");
429  fProfV3Resolution[i]->GetXaxis()->SetBinLabel(8, "<cos(3(#Psi_{TPC} - #Psi_{VZEROC}))>");
430  fProfV3Resolution[i]->GetXaxis()->SetBinLabel(9, "<cos(3(#Psi_{VZERO} - #Psi_{TPC_A}))>");
431  fProfV3Resolution[i]->GetXaxis()->SetBinLabel(10, "<cos(3(#Psi_{VZERO} - #Psi_{TPC_B}))>");
432  fProfV3Resolution[i]->GetXaxis()->SetBinLabel(11, "<cos(3(#Psi_{TPC_A} - #Psi_{TPC_B}))>");
433  fOutput->Add(fProfV3Resolution[i]);
434  fProfV4Resolution[i] = new TProfile(Form("fProfV4Resolution_%i", i), Form("fProfV4Resolution_%i", i), 11, -0.5, 10.5);
435  fProfV4Resolution[i]->GetXaxis()->SetBinLabel(3, "<cos(4(#Psi_{VZEROA} - #Psi_{VZEROC}))>");
436  fProfV4Resolution[i]->GetXaxis()->SetBinLabel(4, "<cos(4(#Psi_{VZEROC} - #Psi_{VZEROA}))>");
437  fProfV4Resolution[i]->GetXaxis()->SetBinLabel(5, "<cos(4(#Psi_{VZEROA} - #Psi_{TPC}))>");
438  fProfV4Resolution[i]->GetXaxis()->SetBinLabel(6, "<cos(4(#Psi_{TPC} - #Psi_{VZEROA}))>");
439  fProfV4Resolution[i]->GetXaxis()->SetBinLabel(7, "<cos(4(#Psi_{VZEROC} - #Psi_{TPC}))>");
440  fProfV4Resolution[i]->GetXaxis()->SetBinLabel(8, "<cos(4(#Psi_{TPC} - #Psi_{VZEROC}))>");
441  fProfV4Resolution[i]->GetXaxis()->SetBinLabel(9, "<cos(4(#Psi_{VZERO} - #Psi_{TPC_A}))>");
442  fProfV4Resolution[i]->GetXaxis()->SetBinLabel(10, "<cos(4(#Psi_{VZERO} - #Psi_{TPC_B}))>");
443  fProfV4Resolution[i]->GetXaxis()->SetBinLabel(11, "<cos(4(#Psi_{TPC_A} - #Psi_{TPC_B}))>");
444  fOutput->Add(fProfV4Resolution[i]);
445  fProfV5Resolution[i] = new TProfile(Form("fProfV5Resolution_%i", i), Form("fProfV5Resolution_%i", i), 11, -0.5, 10.5);
446  fProfV5Resolution[i]->GetXaxis()->SetBinLabel(3, "<cos(5(#Psi_{VZEROA} - #Psi_{VZEROC}))>");
447  fProfV5Resolution[i]->GetXaxis()->SetBinLabel(4, "<cos(5(#Psi_{VZEROC} - #Psi_{VZEROA}))>");
448  fProfV5Resolution[i]->GetXaxis()->SetBinLabel(5, "<cos(5(#Psi_{VZEROA} - #Psi_{TPC}))>");
449  fProfV5Resolution[i]->GetXaxis()->SetBinLabel(6, "<cos(5(#Psi_{TPC} - #Psi_{VZEROA}))>");
450  fProfV5Resolution[i]->GetXaxis()->SetBinLabel(7, "<cos(5(#Psi_{VZEROC} - #Psi_{TPC}))>");
451  fProfV5Resolution[i]->GetXaxis()->SetBinLabel(8, "<cos(5(#Psi_{TPC} - #Psi_{VZEROC}))>");
452  fProfV5Resolution[i]->GetXaxis()->SetBinLabel(9, "<cos(5(#Psi_{VZERO} - #Psi_{TPC_A}))>");
453  fProfV5Resolution[i]->GetXaxis()->SetBinLabel(10, "<cos(5(#Psi_{VZERO} - #Psi_{TPC_B}))>");
454  fProfV5Resolution[i]->GetXaxis()->SetBinLabel(11, "<cos(5(#Psi_{TPC_A} - #Psi_{TPC_B}))>");
455  fOutput->Add(fProfV5Resolution[i]);
456  }
457  }
458 
459  // create histo's used for general QA
460  if (makeQAhistos) {
461  //fHistTPCdEdX = new TH2F("TPCdEdX", "TPCdEdX", 2000, 0.0, 100.0, 500, 0, 500);
462  fHistITSsignal = new TH2F("ITSsignal", "ITSsignal", 2000, 0.0, 100.0, 500, 0, 500);
463  // fHistTOFsignal = new TH2F("TOFsignal", "TOFsignal", 2000, 0.0, 100.0, 500, 0, 500);
464  fHistJetPhi = new TH1F("fHistJetPhi", "Jet #phi Distribution", 72, 0., 1.0*TMath::Pi());
465  fHistTrackPhi = new TH1F("fHistTrackPhi", "Track #phi Distribution", 72, 0., 2.0*TMath::Pi());
466  fHistRhovsCent = new TH2F("RhovsCent", "RhovsCent", 100, 0.0, 100.0, 500, 0, 500);
467  fHistTrackPtallcent = new TH1F("fHistTrackPtallcent", "p_{T} distribution", 1000, 0.0, 100.0);
468  fHistJetEtaPhi = new TH2F("fHistJetEtaPhi", "Jet #eta-#phi",512,-1.8,1.8,512,-3.2,3.2);
469  fHistJetHEtaPhi = new TH2F("fHistJetHEtaPhi", "Jet-Hadron #Delta#eta-#Delta#phi", 72, -1.8, 1.8, 72, -1.6, 4.8);
470  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
471  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
472 
473  fHistClusEtaPhiEnergy = new TH3F("fHistClusEtaPhiEnergy", "Cluster eta-phi-energy", 900, -1.8, 1.8, 720, -3.2, 3.2, 20, 0, 20);
474  fHistClusEnergy = new TH1F("fHistClusEnergy", "Cluster Energy distribution", 200, 0, 20);
475  fHistClusofJetEnergy = new TH1F("fHistClusofJetEnergy", "Cluster of Jet Energy distribution", 200, 0, 20);
476  fHistJetNClusterConstit = new TH1F("fHistJetNClusterConstit", "Number of Cluster Constituents in Jet", 25, 0, 25);
477  fHistJetNTrackConstit = new TH1F("fHistJetNTrackConstit", "Number of Track Constituents in Jet", 25, 0, 25);
478  fHistJetNConstit = new TH1F("fHistJetNConstit", "Number of Total Constituents in Jet", 25, 0, 25);
479 
480  // add to output list
481  //fOutput->Add(fHistTPCdEdX);
482  fOutput->Add(fHistITSsignal);
483  //fOutput->Add(fHistTOFsignal);
484  fOutput->Add(fHistJetPhi);
485  fOutput->Add(fHistTrackPhi);
486  //fOutput->Add(fHistTrackEtaPhi);
488  fOutput->Add(fHistJetEtaPhi);
489  fOutput->Add(fHistJetHEtaPhi);
490  fOutput->Add(fHistSEphieta);
491  fOutput->Add(fHistMEphieta);
493  fOutput->Add(fHistClusEnergy);
498 
499  //for(Int_t ipta=0; ipta<7; ++ipta){
500  // for(Int_t ilab=0; ilab<4; ++ilab){
501  // snprintf(name, nchar, "fHistTrackEtaPhi_%i_%i", ilab,ipta);
502  // fHistTrackEtaPhi[ilab][ipta] = new TH2F(name,name,400,-1,1,640,0.0,2.*TMath::Pi());
503  // fOutput->Add(fHistTrackEtaPhi[ilab][ipta]);
504  // } // end of lab loop
505  //} // end of pt-associated loop
506 
507  for (Int_t i = 0;i<6;++i){
508  name1 = TString(Form("EP0_%i",i));
509  title1 = TString(Form("EP V0 cent bin %i",i));
510  fHistEP0[i] = new TH1F(name1,title1,144,-TMath::Pi(),TMath::Pi());
511  fOutput->Add(fHistEP0[i]);
512 
513  name1 = TString(Form("EP0A_%i",i));
514  title1 = TString(Form("EP V0A cent bin %i",i));
515  fHistEP0A[i] = new TH1F(name1,title1,144,-TMath::Pi(),TMath::Pi());
516  fOutput->Add(fHistEP0A[i]);
517 
518  name1 = TString(Form("EP0C_%i",i));
519  title1 = TString(Form("EP V0C cent bin %i",i));
520  fHistEP0C[i] = new TH1F(name1,title1,144,-TMath::Pi(),TMath::Pi());
521  fOutput->Add(fHistEP0C[i]);
522 
523  name1 = TString(Form("EPAvsC_%i",i));
524  title1 = TString(Form("EP V0A vs V0C cent bin %i",i));
525  fHistEPAvsC[i] = new TH2F(name1,title1,144,-TMath::Pi(),TMath::Pi(),144,-TMath::Pi(),TMath::Pi());
526  fOutput->Add(fHistEPAvsC[i]);
527 
528  name1 = TString(Form("JetPtvsTrackPt_%i",i));
529  title1 = TString(Form("Jet p_{T} vs Leading Track p_{T} cent bin %i",i));
530  fHistJetPtvsTrackPt[i] = new TH2F(name1,title1,250,-50,200,250,0,50);
531  fOutput->Add(fHistJetPtvsTrackPt[i]);
532 
533  name1 = TString(Form("TrackPt_%i",i));
534  title1 = TString(Form("Track p_{T} cent bin %i",i));
535  fHistTrackPt[i] = new TH1F(name1,title1,1000,0,100); // up to 200?
536  fOutput->Add(fHistTrackPt[i]);
537 
538  name1 = TString(Form("JetPtcorrGLrho_%i",i));
539  title1 = TString(Form("Jet p_{T} corrected with Global #rho cent bin %i",i));
540  fHistJetPtcorrGlRho[i] = new TH1F(name1,title1,300,-100,200); // up to 200?
541  fOutput->Add(fHistJetPtcorrGlRho[i]);
542 
543  name1 = TString(Form("JetPtvsdEP_%i",i));
544  title1 = TString(Form("Jet p_{T} vs #DeltaEP cent bin %i",i));
545  fHistJetPtvsdEP[i] = new TH2F(name1,title1,250,-50,200,288,-2*TMath::Pi(),2*TMath::Pi());
546  fOutput->Add(fHistJetPtvsdEP[i]);
547 
548  name1 = TString(Form("RhovsdEP_%i",i));
549  title1 = TString(Form("#rho vs #DeltaEP cent bin %i",i));
550  fHistRhovsdEP[i] = new TH2F(name1,title1,500,0,500,288,-2*TMath::Pi(),2*TMath::Pi());
551  fOutput->Add(fHistRhovsdEP[i]);
552 
553  name1 = TString(Form("JetEtaPhiPt_%i",i));
554  title1 = TString(Form("Jet #eta-#phi p_{T} cent bin %i",i));
555  fHistJetEtaPhiPt[i] = new TH3F(name1,title1,250,-50,200,100,-1,1,64,-3.2,3.2);
556  fOutput->Add(fHistJetEtaPhiPt[i]);
557 
558  name1 = TString(Form("JetPtArea_%i",i));
559  title1 = TString(Form("Jet p_{T} Area cent bin %i",i));
560  fHistJetPtArea[i] = new TH2F(name1,title1,250,-50,200,100,0,1);
561  fOutput->Add(fHistJetPtArea[i]);
562 
563  snprintf(name, nchar, "fHistAreavsRawPt_%i",i);
564  fHistAreavsRawPt[i] = new TH2F(name,name,100,0,1,200,0,200);
565  fOutput->Add(fHistAreavsRawPt[i]);
566  } // loop over centrality
567 
568  } // QA histo switch
569 
570  if (makeBIAShistos) {
571  fHistJetHaddPhiBias = new TH1F("fHistJetHaddPhiBias","Jet-Hadron #Delta#varphi with bias",128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
572  fHistJetHaddPhiINBias = new TH1F("fHistJetHaddPhiINBias","Jet-Hadron #Delta#varphi IN PLANE with bias", 128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
573  fHistJetHaddPhiOUTBias = new TH1F("fHistJetHaddPhiOUTBias","Jet-Hadron #Delta#varphi OUT PLANE with bias",128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
574  fHistJetHaddPhiMIDBias = new TH1F("fHistJetHaddPhiMIDBias","Jet-Hadron #Delta#varphi MIDDLE of PLANE with bias",128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
575 
576  // add to output list
581 
582  for (Int_t i = 0;i<6;++i){
583  name1 = TString(Form("JetPtvsdEPBias_%i",i));
584  title1 = TString(Form("Bias Jet p_{T} vs #DeltaEP cent bin %i",i));
585  fHistJetPtvsdEPBias[i] = new TH2F(name1,title1,250,-50,200,288,-2*TMath::Pi(),2*TMath::Pi());
586  fOutput->Add(fHistJetPtvsdEPBias[i]);
587 
588  name1 = TString(Form("JetEtaPhiPtBias_%i",i));
589  title1 = TString(Form("Jet #eta-#phi p_{T} Bias cent bin %i",i));
590  fHistJetEtaPhiPtBias[i] = new TH3F(name1,title1,250,-50,200,100,-1,1,64,-3.2,3.2);
591  fOutput->Add(fHistJetEtaPhiPtBias[i]);
592 
593  name1 = TString(Form("JetPtAreaBias_%i",i));
594  title1 = TString(Form("Jet p_{T} Area Bias cent bin %i",i));
595  fHistJetPtAreaBias[i] = new TH2F(name1,title1,250,-50,200,100,0,1);
596  fOutput->Add(fHistJetPtAreaBias[i]);
597  } // end of centrality loop
598  } // bias histos
599 
600  if (makeoldJEThadhistos) {
601  for(Int_t icent = 0; icent<6; ++icent){
602  snprintf(name, nchar, "fHistJetPtTT_%i",icent);
603  fHistJetPtTT[icent] = new TH1F(name,name,200,0,200);
604  fOutput->Add(fHistJetPtTT[icent]);
605 
606  snprintf(name, nchar, "fHistJetPt_%i",icent);
607  fHistJetPt[icent] = new TH1F(name,name,200,0,200);
608  fOutput->Add(fHistJetPt[icent]);
609 
610  snprintf(name, nchar, "fHistJetPtBias_%i",icent);
611  fHistJetPtBias[icent] = new TH1F(name,name,200,0,200);
612  fOutput->Add(fHistJetPtBias[icent]);
613 
614  for(Int_t iptjet = 0; iptjet<5; ++iptjet){
615  for(Int_t ieta = 0; ieta<3; ++ieta){
616  snprintf(name, nchar, "fHistJetH_%i_%i_%i",icent,iptjet,ieta);
617  fHistJetH[icent][iptjet][ieta]=new TH2F(name,name,64,-0.5*TMath::Pi(),1.5*TMath::Pi(),300,0,30);
618  fOutput->Add(fHistJetH[icent][iptjet][ieta]);
619 
620  snprintf(name, nchar, "fHistJetHBias_%i_%i_%i",icent,iptjet,ieta);
621  fHistJetHBias[icent][iptjet][ieta]=new TH2F(name,name,64,-0.5*TMath::Pi(),1.5*TMath::Pi(),300,0,30);
622  fOutput->Add(fHistJetHBias[icent][iptjet][ieta]);
623 
624  snprintf(name, nchar, "fHistJetHTT_%i_%i_%i",icent,iptjet,ieta);
625  fHistJetHTT[icent][iptjet][ieta]=new TH2F(name,name,64,-0.5*TMath::Pi(),1.5*TMath::Pi(),300,0,30);
626  fOutput->Add(fHistJetHTT[icent][iptjet][ieta]);
627  } // end of eta loop
628  } // end of pt-jet loop
629  } // end of centrality loop
630  } // make JetHadhisto old
631 
632  if (makeextraCORRhistos) {
633  for(Int_t itrackpt=0; itrackpt<9; itrackpt++){
634  snprintf(name, nchar, "fHistJetHadbindPhi_%i",itrackpt);
635  fHistJetHadbindPhi[itrackpt] = new TH1F(name,name,128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
636  fOutput->Add(fHistJetHadbindPhi[itrackpt]);
637 
638  snprintf(name, nchar, "fHistJetHadbindPhiIN_%i",itrackpt);
639  fHistJetHadbindPhiIN[itrackpt] = new TH1F(name,name,128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
640  fOutput->Add(fHistJetHadbindPhiIN[itrackpt]);
641 
642  snprintf(name, nchar, "fHistJetHadbindPhiMID_%i",itrackpt);
643  fHistJetHadbindPhiMID[itrackpt] = new TH1F(name,name,128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
644  fOutput->Add(fHistJetHadbindPhiMID[itrackpt]);
645 
646  snprintf(name, nchar, "fHistJetHadbindPhiOUT_%i",itrackpt);
647  fHistJetHadbindPhiOUT[itrackpt] = new TH1F(name,name,128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
648  fOutput->Add(fHistJetHadbindPhiOUT[itrackpt]);
649  } // end of trackpt bin loop
650 
651  for (Int_t i = 0;i<6;++i){
652  name1 = TString(Form("JetHaddPhiINcent_%i",i));
653  title1 = TString(Form("Jet Hadron #Delta#varphi Distribution IN PLANE cent bin %i",i));
654  fHistJetHaddPhiINcent[i] = new TH1F(name1,title1,128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
656 
657  name1 = TString(Form("JetHaddPhiOUTcent_%i",i));
658  title1 = TString(Form("Jet Hadron #Delta#varphi Distribution OUT PLANE cent bin %i",i));
659  fHistJetHaddPhiOUTcent[i] = new TH1F(name1,title1,128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
661 
662  name1 = TString(Form("JetHaddPhiMIDcent_%i",i));
663  title1 = TString(Form("Jet Hadron #Delta#varphi Distribution MIDDLE of PLANE cent bin %i",i));
664  fHistJetHaddPhiMIDcent[i] = new TH1F(name1,title1,128,-0.5*TMath::Pi(), 1.5*TMath::Pi());
666 
667  name1 = TString(Form("JetPtNcon_%i",i));
668  title1 = TString(Form("Jet p_{T} Ncon cent bin %i",i));
669  fHistJetPtNcon[i] = new TH2F(name1,title1,250,-50,200,100,0,2000);
670  fOutput->Add(fHistJetPtNcon[i]);
671 
672  name1 = TString(Form("JetPtNconBias_%i",i));
673  title1 = TString(Form("Jet p_{T} NconBias cent bin %i",i));
674  fHistJetPtNconBias[i] = new TH2F(name1,title1,250,-50,200,100,0,2000);
675  fOutput->Add(fHistJetPtNconBias[i]);
676 
677  name1 = TString(Form("JetPtNconCh_%i",i));
678  title1 = TString(Form("Jet p_{T} NconCh cent bin %i",i));
679  fHistJetPtNconCh[i] = new TH2F(name1,title1,250,-50,200,100,0,2000);
680  fOutput->Add(fHistJetPtNconCh[i]);
681 
682  name1 = TString(Form("JetPtNconBiasCh_%i",i));
683  title1 = TString(Form("Jet p_{T} NconBiasCh cent bin %i",i));
684  fHistJetPtNconBiasCh[i] = new TH2F(name1,title1,250,-50,200,100,0,2000);
685  fOutput->Add(fHistJetPtNconBiasCh[i]);
686 
687  name1 = TString(Form("JetPtNconEm_%i",i));
688  title1 = TString(Form("Jet p_{T} NconEm cent bin %i",i));
689  fHistJetPtNconEm[i] = new TH2F(name1,title1,250,-50,200,100,0,2000);
690  fOutput->Add(fHistJetPtNconEm[i]);
691 
692  name1 = TString(Form("JetPtNconBiasEm_%i",i));
693  title1 = TString(Form("Jet p_{T} NconBiasEm cent bin %i",i));
694  fHistJetPtNconBiasEm[i] = new TH2F(name1,title1,250,-50,200,100,0,2000);
695  fOutput->Add(fHistJetPtNconBiasEm[i]);
696  } // extra Correlations histos switch
697  }
698 
699  // variable binned pt for THnSparse's
700  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};
701  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};
702 
703  // tracks: 51, jets: 26
704  // number of bins you tell histogram should be (# in array - 1) because the last bin
705  // is the right-most edge of the histogram
706  // i.e. this is for PT and there are 57 numbers (bins) thus we have 56 bins in our histo
707  Int_t nbinsjetPT = sizeof(xlowjetPT)/sizeof(Double_t) - 1;
708  Int_t nbinstrPT = sizeof(xlowtrPT)/sizeof(Double_t) - 1;
709 
710  // set up jet-hadron sparse
711  UInt_t bitcodeMESE = 0; // bit coded, see GetDimParams() below
712  UInt_t bitcodePID = 0; // bit coded, see GetDimParamsPID() below
713  UInt_t bitcodeCorr = 0; // bit coded, see GetDimparamsCorr() below
714  bitcodeMESE = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6; // | 1<<7 | 1<<8 | 1<<9;
715  if(fDoEventMixing) {
716  fhnJH = NewTHnSparseF("fhnJH", bitcodeMESE);
717 
718  if(dovarbinTHnSparse){
719  fhnJH->GetAxis(1)->Set(nbinsjetPT, xlowjetPT);
720  fhnJH->GetAxis(2)->Set(nbinstrPT, xlowtrPT);
721  }
722 
723  fOutput->Add(fhnJH);
724  }
725 
726  bitcodeCorr = 1<<0 | 1<<1 | 1<<2 | 1<<3; // | 1<<4 | 1<<5;
727  fhnCorr = NewTHnSparseFCorr("fhnCorr", bitcodeCorr);
728  if(dovarbinTHnSparse) fhnCorr->GetAxis(1)->Set(nbinsjetPT, xlowjetPT);
729  fOutput->Add(fhnCorr);
730 
731  /*
732  Double_t centralityBins[nCentralityBins+1];
733  for(Int_t ic=0; ic<nCentralityBins+1; ic++){
734  if(ic==nCentralityBins) centralityBins[ic]=500;
735  else centralityBins[ic]=10.0*ic;
736  }
737  */
738 
739  // set up centrality bins for mixed events
740  // for pp we need mult bins for event mixing. Create binning here, to also make a histogram from it
741  Int_t nCentralityBinspp = 8;
742  //Double_t centralityBinspp[nCentralityBinspp+1];
743  Double_t centralityBinspp[9] = {0.0, 4., 9, 15, 25, 35, 55, 100.0, 500.0};
744 
745  // Setup for Pb-Pb collisions
746  Int_t nCentralityBinsPbPb = 100;
747  Double_t mult = 1.0;
748  if(fCentBinSize==1) {
749  nCentralityBinsPbPb = 100;
750  mult = 1.0;
751  } else if(fCentBinSize==2){
752  nCentralityBinsPbPb = 50;
753  mult = 2.0;
754  } else if(fCentBinSize==5){
755  nCentralityBinsPbPb = 20;
756  mult = 5.0;
757  } else if(fCentBinSize==10){
758  nCentralityBinsPbPb = 10;
759  mult = 10.0;
760  }
761 
762  Double_t centralityBinsPbPb[nCentralityBinsPbPb]; // nCentralityBinsPbPb
763  for(Int_t ic=0; ic<nCentralityBinsPbPb; ic++){
764  centralityBinsPbPb[ic]=mult*ic;
765  }
766 
767 /*
768  Int_t nCentralityBinsPbPb = 10; //100;
769  Double_t centralityBinsPbPb[nCentralityBinsPbPb+1];
770  for(Int_t ic=0; ic<nCentralityBinsPbPb; ic++){
771  centralityBinsPbPb[ic]=10.0*ic; //1.0*ic;
772  }
773 */
774 
775  if(fBeam == 0) fHistMult = new TH1F("fHistMult","multiplicity",nCentralityBinspp,centralityBinspp);
776  if(fBeam == 1) fHistMult = new TH1F("fHistMult","multiplicity",nCentralityBinsPbPb,centralityBinsPbPb);
777 // fOutput->Add(fHistMult);
778 
779  // Event Mixing
780  Int_t trackDepth = fMixingTracks;
781  Int_t poolsize = 1000; // Maximum number of events, ignored in the present implemented of AliEventPoolManager
782  Int_t nZvtxBins = 5+1+5;
783  Double_t vertexBins[] = {-10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10};
784  Double_t* zvtxbin = vertexBins;
785  if(fBeam == 0) fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBinspp, centralityBinspp, nZvtxBins, zvtxbin);
786  if(fBeam == 1) fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBinsPbPb, centralityBinsPbPb, nZvtxBins, zvtxbin);
787 
788  // set up event mixing sparse
789  if(fDoEventMixing){
790  bitcodeMESE = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6; // | 1<<7 | 1<<8 | 1<<9;
791  fhnMixedEvents = NewTHnSparseF("fhnMixedEvents", bitcodeMESE);
792 
793  // set up some variable binning of the sparse
794  if(dovarbinTHnSparse){
795  fhnMixedEvents->GetAxis(1)->Set(nbinsjetPT, xlowjetPT);
796  fhnMixedEvents->GetAxis(2)->Set(nbinstrPT, xlowtrPT);
797  }
798 
799  fOutput->Add(fhnMixedEvents);
800  } // end of do-eventmixing
801 
802  // set up PID sparse
803  if(doPID){
804  // ****************************** PID *****************************************************
805  // set up PID handler
806  AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
807  AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
808  if(!inputHandler) {
809  AliFatal("Input handler needed");
810  return;
811  }
812 
813  // PID response object
814  fPIDResponse = inputHandler->GetPIDResponse();
815  if (!fPIDResponse) {
816  AliError("PIDResponse object was not created");
817  return;
818  }
819  // *****************************************************************************************
820 
821  // PID counter
822  fHistPID = new TH1F("fHistPID","PID Counter", 15, 0.5, 25.5);
824  fOutput->Add(fHistPID);
825 
826  if(allpidAXIS) {
827  bitcodePID = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<7 | 1<<8 | 1<<9 |
828  1<<10 | 1<<11 | 1<<12 | 1<<13 | 1<<14 | 1<<15 | 1<<16 | 1<<17 | 1<<18; // | 1<<19 | 1<<20;
829  fhnPID = NewTHnSparseFPID("fhnPID", bitcodePID);
830  } else {
831  bitcodePID = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<7 | 1<<8 | 1<<9 |
832  1<<10 | 1<<11; // | 1<<12 | 1<<13;
833  fhnPID = NewTHnSparseFPID("fhnPID", bitcodePID);
834  }
835 
836  // set some variable binning of sparse
837  if(dovarbinTHnSparse){
838  fhnPID->GetAxis(1)->Set(nbinstrPT, xlowtrPT);
839  fhnPID->GetAxis(8)->Set(nbinsjetPT, xlowjetPT);
840  }
841 
842  fOutput->Add(fhnPID);
843  } // end of do-PID
844 
845  // =========== Switch on Sumw2 for all histos ===========
846  for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
847  TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
848  if (h1){
849  h1->Sumw2();
850  continue;
851  }
852  TH2 *h2 = dynamic_cast<TH2*>(fOutput->At(i));
853  if (h2){
854  h2->Sumw2();
855  continue;
856  }
857  TH3 *h3 = dynamic_cast<TH3*>(fOutput->At(i));
858  if (h3){
859  h3->Sumw2();
860  continue;
861  }
862  THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
863  if(hn)hn->Sumw2();
864  }
865 
866  // Added March14, 2016 - used for new framework (track container)
868  fTracksCont->SetClassName("AliVTrack");
869 
870  PostData(1, fOutput);
871 }
872 
873 //_________________________________________________________________________
875 {
877 
878  // Added March14, 2016 - used for new framework, to filled ClonesArray with container contents
879  fTracksFromContainer = new TClonesArray("AliVParticle");
880  fTracksFromContainer->SetName("TracksFromContainer");
881  fTracksFromContainer->SetOwner(kTRUE);
882 
883  if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
884  if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
886 }
887 
888 //_____________________________________________________________________________
890 {
891  // determine the run number to see if the track and jet cuts should be refreshed for semi-good TPC runs
893  return kTRUE;
894 }
895 
896 //_________________________________________________________________________
898 { // Main loop called for each event
899  // TEST TEST TEST TEST for OBJECTS!
900 
901  fHistEventQA->Fill(1); // All Events that get entered
902 
903  // check and fill a Event Selection QA histogram for different trigger selections
904  UInt_t trig = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
905 
906  // fill Event Trigger QA
908 
909  // see if we are running on PbPb and try and get LocalRho object
910  if(GetBeamType() == 1) {
911  if(!fLocalRho){
912  AliError(Form("Couldn't get fLocalRho object, try to get it from Event based on name\n"));
914  if(!fLocalRho) return kTRUE;
915  }
916  } // check for LocalRho object if PbPb data
917 
918  if(!fTracks){
919  AliError(Form("No fTracks object!!\n"));
920  return kTRUE;
921  }
922  if(!fJets){
923  AliError(Form("No fJets object!!\n"));
924  //return kTRUE;
925  }
926 
927  fHistEventQA->Fill(2); // events after object check
928 
929  // what kind of event do we have: AOD or ESD?
930  Bool_t useAOD;
931  if (dynamic_cast<AliAODEvent*>(InputEvent())) useAOD = kTRUE;
932  else useAOD = kFALSE;
933 
934  // if we have ESD event, set up ESD object
935  if(!useAOD){
936  fESD = dynamic_cast<AliESDEvent*>(InputEvent());
937  if (!fESD) {
938  AliError(Form("ERROR: fESD not available\n"));
939  return kTRUE;
940  }
941  }
942 
943  // if we have AOD event, set up AOD object
944  if(useAOD){
945  fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
946  if(!fAOD) {
947  AliError(Form("ERROR: fAOD not available\n"));
948  return kTRUE;
949  }
950  }
951 
952  fHistEventQA->Fill(3); // events after Aod/esd check
953 
954  // get centrality
955  Int_t centbin = GetCentBin(fCent);
956 
957  // to limit filling unused entries in sparse, only fill for certain centrality ranges
958  // ranges can be different than functional cent bin setter
959  Int_t cbin = -1;
960  if (fCent>=0 && fCent<10) cbin = 1;
961  else if (fCent>=10 && fCent<20) cbin = 2;
962  else if (fCent>=20 && fCent<30) cbin = 3;
963  else if (fCent>=30 && fCent<50) cbin = 4;
964  else if (fCent>=50 && fCent<90) cbin = 5;
965  else cbin = -99;
966 
968 /*
969  if(fCent <= fCentralityClasses->At(0) || fCent >= fCentralityClasses->At(fCentralityClasses->GetSize()-1) || TMath::Abs(fCent-InputEvent()->GetCentrality()->GetCentralityPercentile("TRK")) > 5.) return kFALSE;
970  // determine centrality class
971  fInCentralitySelection = -1;
972  for(Int_t i(0); i < fCentralityClasses->GetSize()-1; i++) {
973  if(fCent >= fCentralityClasses->At(i) && fCent <= fCentralityClasses->At(1+i)) {
974  fInCentralitySelection = i;
975  break;
976  }
977  }
978  if(fInCentralitySelection<0) return kFALSE; // should be null op
980 */
981 
982  // if we are on PbPb data do cut on centrality > 90%, else by default DON'T
983  if (GetBeamType() == 1) {
984  // apply cut to event on Centrality > 90%
985  if(fCent>90) return kTRUE;
986  }
987 
988  // BEAM TYPE enumerator: kNA = -1, kpp = 0, kAA = 1, kpA = 2
989  // for pp analyses we will just use the first centrality bin
990  if(GetBeamType() == 0) if (centbin == -1) centbin = 0;
991  if(GetBeamType() == 1) if (centbin == -1) return kTRUE;
992 
993  fHistEventQA->Fill(4); // events after centrality check
994 
995  // get vertex information
996  Double_t fvertex[3]={0,0,0};
997  InputEvent()->GetPrimaryVertex()->GetXYZ(fvertex);
998  Double_t zVtx=fvertex[2];
999 
1000  // get z-vertex bin
1001  //Int_t zVbin = GetzVertexBin(zVtx);
1002 
1003  // apply zVtx cut
1004  if(fabs(zVtx)>10.0) return kTRUE;
1005 
1006  fHistEventQA->Fill(5); // events after zvertex check
1007 
1008  // create pointer to list of input event
1009  TList *list = InputEvent()->GetList();
1010  if(!list) {
1011  AliError(Form("ERROR: list not attached\n"));
1012  return kTRUE;
1013  }
1014 
1015  fHistEventQA->Fill(6); // events after list check
1016 
1017 //test
1018  GetTriggerList();
1019 
1020  // initialize TClonesArray pointers to jets and tracks
1021  TClonesArray *jets = 0;
1022  TClonesArray *tracks = 0;
1023  TClonesArray *tracksME = 0;
1024  TClonesArray *clusters = 0;
1025 
1026  // method for filling collection output array of 'saved' tracks - Added March14, 2016 for new framework
1027  Int_t tacc = 0;
1028  fTracksFromContainer->Clear();
1029  if(fTracksCont) {
1031  AliVTrack *track = dynamic_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
1032  while(track) {
1033  // add container tracks to clones array
1034  (*fTracksFromContainer)[tacc] = track;
1035  ++tacc;
1036 
1037  // get next track
1038  track = dynamic_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
1039  } // accepted track loop
1040  }
1041 
1042  // get Tracks object
1043  tracks = fTracksFromContainer; // added March14, 2016
1044  if (!tracks) {
1045  AliError(Form("Pointer to tracks: %s == 0", fTracksCont->GetName()));
1046  return kTRUE;
1047  } // verify existence of tracks
1048 
1049  // get ME Tracks object
1050  tracksME = fTracksFromContainer; // Added March14, 2016
1051  if (!tracksME) {
1052  AliError(Form("Pointer to ME tracks: %s == 0", fTracksCont->GetName()));
1053  return kTRUE;
1054  } // verify existence of tracks
1055 
1056  // get Jets object
1057  Int_t Njets = 0;
1058  if(fJets) {
1059  jets = dynamic_cast<TClonesArray*>(list->FindObject(fJets));
1060  if(!jets){
1061  AliError(Form("Pointer to jets: %s == 0", fJets->GetName()));
1062  //return kTRUE;
1063  } // verify existence of jets
1064  Njets = jets->GetEntries();
1065  }
1066 
1067  fHistEventQA->Fill(7); // events after track/jet pointer check
1068 
1069  // get number of jets and tracks
1070 // Int_t Njets;
1071 // if(!fJets) { Njets = 0; }
1072 // else Njets = jets->GetEntries();
1073  //const Int_t Njets = jets->GetEntries();
1074  const Int_t Ntracks = tracks->GetEntries();
1077 
1078  fHistEventQA->Fill(8); // events after #track and jets < 1 check
1079 
1080  if (makeQAhistos) fHistMult->Fill(Ntracks); // fill multiplicity distribution
1081 
1082  // get cluster collection (mainly QA at this point)
1083  clusters = dynamic_cast<TClonesArray*>(list->FindObject(fCaloClustersName));
1084  if (!clusters) {
1085  AliError(Form("Pointer to clusters: %s == 0", fCaloClustersName.Data()));
1086  return kTRUE;
1087  } // verify existence of clusters
1088 
1089  // get clusters and loop over
1090  const Int_t Nclusters = clusters->GetEntries();
1091  for (Int_t iclus = 0; iclus < Nclusters; iclus++){
1092  AliVCluster* cluster = static_cast<AliVCluster*>(clusters->At(iclus));
1093  if(!cluster){
1094  AliError(Form("Couldn't get AliVCluster %d\n", iclus));
1095  continue;
1096  }
1097 
1098  // get some info on cluster and fill histo
1099  TLorentzVector nPart;
1100  cluster->GetMomentum(nPart, fvertex);
1101  if(makeQAhistos) fHistClusEtaPhiEnergy->Fill(nPart.Eta(), nPart.Phi(), nPart.E());
1102  if(makeQAhistos) fHistClusEnergy->Fill(nPart.E());
1103  }
1104 
1105  // initialize track parameters
1106  Int_t iTT=-1;
1107  Double_t ptmax=-10;
1108  Int_t NtrackAcc = 0;
1109 
1110  fVevent = dynamic_cast<AliVEvent*>(InputEvent());
1111  if (!fVevent) {
1112  printf("ERROR: fVEvent not available\n");
1113  return kTRUE;
1114  }
1115 
1116  // fill event mixing QA
1117  if(trig & AliVEvent::kEMCEGA) fHistCentZvertGA->Fill(fCent, zVtx);
1118  if(trig & AliVEvent::kEMCEJE) fHistCentZvertJE->Fill(fCent, zVtx);
1119  if(trig & AliVEvent::kMB) fHistCentZvertMB->Fill(fCent, zVtx);
1120  if(trig & AliVEvent::kAny) fHistCentZvertAny->Fill(fCent, zVtx);
1121 
1122  // loop over tracks - to get hardest track (highest pt)
1123  for (Int_t iTracks = 0; iTracks < Ntracks; iTracks++){
1124  AliVTrack *track = dynamic_cast<AliVTrack*>(tracks->At(iTracks));
1125  if(!track){
1126  AliError(Form("Couldn't get AliVTrack %d\n", iTracks));
1127  continue;
1128  }
1129 
1130  // track cuts
1131  if(TMath::Abs(track->Eta())>0.9) continue;
1132  if(track->Pt()<0.15) continue;
1133  //iCount++;
1134  NtrackAcc++;
1135 
1136  if(track->Pt()>ptmax){
1137  ptmax=track->Pt(); // max pT track
1138  iTT=iTracks; // trigger tracks
1139  } // check if Pt>maxpt
1140 
1141  if (makeQAhistos) fHistTrackPhi->Fill(track->Phi());
1142  if (makeQAhistos) fHistTrackPt[centbin]->Fill(track->Pt());
1143  if (makeQAhistos) fHistTrackPtallcent->Fill(track->Pt());
1144  } // end of loop over tracks
1145 
1146  if(doEventPlaneRes){
1147  // cut on event selection before calculating reaction plane and filling histo's
1148  if ((trig && fTriggerEventType)) {
1149  // cache the leading jet within acceptance
1151 
1152  // set storage vectors for 2nd and 3rd order event plane values for different subevents
1153  Double_t vzero[2][2];
1154  /* for the combined vzero event plane
1155  * [0] psi2 [1] psi3
1156  */
1157  Double_t vzeroComb[2];
1158  // [0] psi2 [1] psi3
1159  Double_t tpc[2];
1160  // evaluate the actual event planes
1161  // grab the actual data
1162  CalculateEventPlaneVZERO(vzero);
1165  CalculateEventPlaneResolution(vzero, vzeroComb, tpc);
1166  }
1167  }
1168 
1169  // get rho from event and fill relative histo's
1171  fRhoVal = fRho->GetVal();
1172 
1173  if (makeQAhistos) {
1174  fHistRhovsdEP[centbin]->Fill(fRhoVal,fEPV0); // Global Rho vs delta Event Plane angle
1175  fHistRhovsCent->Fill(fCent,fRhoVal); // Global Rho vs Centrality
1176  fHistEP0[centbin]->Fill(fEPV0);
1177  fHistEP0A[centbin]->Fill(fEPV0A);
1178  fHistEP0C[centbin]->Fill(fEPV0C);
1179  fHistEPAvsC[centbin]->Fill(fEPV0A,fEPV0C);
1180  }
1181 
1182  // initialize jet parameters
1183  Int_t ijethi=-1;
1184  Double_t highestjetpt=0.0;
1185  Int_t passedTTcut=0;
1186  Int_t NjetAcc = 0;
1187  Double_t leadhadronPT = 0.0;
1188  Double_t maxclusterpt = 0.0;
1189  Double_t maxtrackpt = 0.0;
1190 
1191  // loop over jets in an event - to find highest jet pT and apply some cuts
1192  for (Int_t ijet = 0; ijet < Njets; ijet++){
1193  // get our jets
1194  AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(ijet));
1195  if (!jet) continue;
1196 
1197  // apply jet cuts
1198  if ((jet->Phi()<fPhimin)||(jet->Phi()>fPhimax)) continue;
1199  if ((jet->Eta()<fEtamin)||(jet->Eta()>fEtamax)) continue;
1200  if (makeQAhistos) fHistAreavsRawPt[centbin]->Fill(jet->Pt(),jet->Area());
1201  if(!AcceptMyJet(jet)) continue;
1202 
1203  NjetAcc++; // # of accepted jets
1204 
1205  // if FlavourJetAnalysis, get desired FlavTag and check against Jet
1206  if(doFlavourJetAnalysis) { if(!AcceptFlavourJet(jet, fJetFlavTag)) continue;}
1207 
1208  // use this to get total # of jets passing cuts in events!!!!!!!!
1209  if (makeQAhistos) fHistJetPhi->Fill(jet->Phi()); // Jet Phi histogram (filled)
1210 
1211  // get highest Pt jet in event
1212  if(highestjetpt<jet->Pt()){
1213  ijethi=ijet;
1214  highestjetpt=jet->Pt();
1215  }
1216  } // end of looping over jets
1217 
1218  // accepted jets
1219  fHistNjetvsCent->Fill(fCent,NjetAcc);
1220  Int_t NJETAcc = 0;
1221  fHistEventQA->Fill(9); // events after track/jet loop to get highest pt
1222 
1223  //cout<<"Event #: "<<event<<endl;
1224 
1225  // loop over jets in event and make appropriate cuts
1226  for (Int_t ijet = 0; ijet < Njets; ++ijet) {
1227  AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(ijet));
1228  if (!jet) continue;
1229 
1230  // check our jet is in our selected trigger event
1231  if (!(trig && fTriggerEventType)) continue;
1232 
1233  // (should probably be higher..., but makes a cut on jet pT)
1234  if (jet->Pt()<0.1) continue;
1235  // do we accept jet? apply jet cuts
1236  if (!AcceptMyJet(jet)) continue;
1237 
1238  // if FlavourJetAnalysis, get desired FlavTag and check against Jet
1239  if(doFlavourJetAnalysis) { if(!AcceptFlavourJet(jet, fJetFlavTag)) continue;}
1240 
1241  fHistEventQA->Fill(10); // accepted jets
1242 
1243  // check on lead jet
1244  Double_t leadjet=0;
1245  if (ijet==ijethi) leadjet=1;
1246 
1247  // check on leading hadron pt
1248  if (ijet==ijethi) leadhadronPT = GetLeadingHadronPt(jet);
1249  if (ijet==ijethi) maxclusterpt = jet->MaxClusterPt();
1250  if (ijet==ijethi) maxtrackpt = jet->MaxTrackPt();
1251 
1252  // initialize and calculate various parameters: pt, eta, phi, rho, etc...
1253  Double_t jetphi = jet->Phi(); // phi of jet
1254  NJETAcc++; // # accepted jets
1255  Double_t jeteta = jet->Eta(); // ETA of jet
1256  Double_t jetPt = -500;
1257  Double_t jetPtGlobal = -500;
1258  Double_t jetPtLocal = -500; // initialize corr jet pt
1259  if(GetBeamType() == 1) {
1260  fLocalRhoVal = fLocalRho->GetLocalVal(jetphi, fJetRad); //GetJetRadius(0)); // get local rho value
1261  jetPtLocal = jet->Pt()-jet->Area()*fLocalRhoVal; // corrected pT of jet using Rho modulated for V2 and V3
1262  }
1263  jetPt = jet->Pt();
1264  jetPtGlobal = jet->Pt()-jet->Area()*fRhoVal; // corrected pT of jet from rho value
1265  Double_t dEP = -500; // initialize angle between jet and event plane
1266  dEP = RelativeEPJET(jetphi,fEPV0); // angle betweeen jet and event plane
1267 
1268  // make histo's
1269  if(makeQAhistos) fHistJetPtvsTrackPt[centbin]->Fill(jetPt,jet->MaxTrackPt());
1270  if(makeQAhistos) fHistJetPtcorrGlRho[centbin]->Fill(jetPtGlobal);
1271  if(makeQAhistos) fHistJetPtvsdEP[centbin]->Fill(jetPt, dEP);
1272  if(makeQAhistos) fHistJetEtaPhiPt[centbin]->Fill(jetPt,jet->Eta(),jet->Phi());
1273  if(makeQAhistos) fHistJetPtArea[centbin]->Fill(jetPt,jet->Area());
1274  if(makeQAhistos) fHistJetEtaPhi->Fill(jet->Eta(),jet->Phi()); // fill jet eta-phi distribution histo
1275  if(makeextraCORRhistos) fHistJetPtNcon[centbin]->Fill(jetPt,jet->GetNumberOfConstituents());
1276  if(makeextraCORRhistos) fHistJetPtNconCh[centbin]->Fill(jetPt,jet->GetNumberOfTracks());
1277  if(makeextraCORRhistos) fHistJetPtNconEm[centbin]->Fill(jetPt,jet->GetNumberOfClusters());
1278  if(makeoldJEThadhistos) fHistJetPt[centbin]->Fill(jet->Pt()); // fill #jets vs pT histo
1279  //fHistDeltaPtvsArea->Fill(jetPt,jet->Area());
1280 
1281  // make histo's with BIAS applied
1282  if (jet->MaxTrackPt()>fTrkBias){
1283  if(makeBIAShistos) fHistJetPtvsdEPBias[centbin]->Fill(jetPt, dEP);
1284  if(makeBIAShistos) fHistJetEtaPhiPtBias[centbin]->Fill(jetPt,jet->Eta(),jet->Phi());
1285  if(makeextraCORRhistos) fHistJetPtAreaBias[centbin]->Fill(jetPt,jet->Area());
1286  if(makeextraCORRhistos) fHistJetPtNconBias[centbin]->Fill(jetPt,jet->GetNumberOfConstituents());
1287  if(makeextraCORRhistos) fHistJetPtNconBiasCh[centbin]->Fill(jetPt,jet->GetNumberOfTracks());
1288  if(makeextraCORRhistos) fHistJetPtNconBiasEm[centbin]->Fill(jetPt,jet->GetNumberOfClusters());
1289  }
1290 
1291  //if(leadjet && centbin==0){
1292  // if(makeextraCORRhistos) fHistJetPt[centbin+1]->Fill(jet->Pt());
1293  //}
1294  if ((jet->MaxTrackPt()>fTrkBias) || (jet->MaxClusterPt()>fClusBias)){
1295  if (makeoldJEThadhistos){
1296  fHistJetPtBias[centbin]->Fill(jet->Pt());
1297  //if(leadjet && centbin==0) fHistJetPtBias[centbin+1]->Fill(jet->Pt());
1298  }
1299  } // end of MaxTrackPt>ftrkBias or maxclusterPt>fclusBias
1300 
1301  // do we have trigger tracks
1302  if(iTT>0){
1303  AliVTrack* TT = static_cast<AliVTrack*>(tracks->At(iTT));
1304  if(TMath::Abs(jet->Phi()-TT->Phi()-TMath::Pi())<0.6) passedTTcut=1;
1305  else passedTTcut=0;
1306  } // end of check on iTT > 0
1307 
1308  if(passedTTcut) {
1309  if (makeoldJEThadhistos) fHistJetPtTT[centbin]->Fill(jet->Pt());
1310  }
1311 
1312  // cut on HIGHEST jet pt in event (15 GeV default)
1313  //if (highestjetpt>fJetPtcut) {
1314  if (jet->Pt() > fJetPtcut) {
1315  fHistEventQA->Fill(11); // jets meeting pt threshold
1316 
1317  // does our max track or cluster pass the bias?
1318  if ((jet->MaxTrackPt()>fTrkBias) || (jet->MaxClusterPt()>fClusBias)){
1319  // set up and fill Jet-Hadron trigger jets THnSparse
1320  if(fcorrJetPt) {
1321  Double_t CorrEntries[4] = {fCent, jetPtLocal, dEP, zVtx};
1322  if(fReduceStatsCent > 0) {
1323  if(cbin == fReduceStatsCent) fhnCorr->Fill(CorrEntries); // fill Sparse Histo with trigger Jets entries
1324  } else fhnCorr->Fill(CorrEntries); // fill Sparse Histo with trigger Jets entries
1325  } else { // don't correct jet pt
1326  Double_t CorrEntries[4] = {fCent, jet->Pt(), dEP, zVtx};
1327  if(fReduceStatsCent > 0) {
1328  if(cbin == fReduceStatsCent) fhnCorr->Fill(CorrEntries); // fill Sparse Histo with trigger Jets entries
1329  } else fhnCorr->Fill(CorrEntries); // fill Sparse Histo with trigger Jets entries
1330  }
1331 
1332  if(GetBeamType() == 1) fHistLocalRhoJetpt->Fill(jetPtLocal);
1336 
1337  // get clusters and loop over
1338  for (Int_t iclus = 0; iclus < Nclusters; iclus++){
1339  AliVCluster* clusterofJet = static_cast<AliVCluster*>(clusters->At(iclus));
1340  if(!clusterofJet){
1341  AliError(Form("Couldn't get AliVCluster %d\n", iclus));
1342  continue;
1343  }
1344 
1345  // get some info on cluster and fill histo
1346  TLorentzVector nPart;
1347  clusterofJet->GetMomentum(nPart, fvertex);
1348  //if(makeQAhistos) fHistClusEtaPhiEnergy->Fill(nPart.Eta(), nPart.Phi(), nPart.E());
1349  if(makeQAhistos) fHistClusofJetEnergy->Fill(nPart.E());
1350 
1351  } // loop over clusters for each jet
1352  } // check on max track and cluster pt/Et
1353 
1354  // loop over all track for an event containing a jet with a pT>fJetPtCut (15)GeV
1355  for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
1356  AliVTrack *track = dynamic_cast<AliVTrack*>(tracks->At(iTracks));
1357  if (!track) {
1358  AliError(Form("Couldn't get AliVtrack %d\n", iTracks));
1359  continue;
1360  }
1361 
1362  // apply track cuts
1363  if(TMath::Abs(track->Eta())>fTrkEta) continue;
1364  if (track->Pt()<0.15) continue;
1365 
1366  fHistEventQA->Fill(12); // accepted tracks in events from trigger jets
1367 
1368  // calculate and get some track parameters
1369  Double_t trCharge = -99;
1370  trCharge = track->Charge();
1371  Double_t tracketa=track->Eta(); // eta of track
1372  Double_t deta=tracketa-jeteta; // dETA between track and jet
1373  //Double_t dR=sqrt(deta*deta+dphijh*dphijh); // difference of R between jet and hadron track
1374 
1375  // calculate single particle tracking efficiency for correlations
1376  Double_t trefficiency = -999;
1377  trefficiency = EffCorrection(track->Eta(), track->Pt(), fDoEffCorr);
1378 
1379  Int_t ieta = -1; // initialize deta bin
1380  Int_t iptjet = -1; // initialize jet pT bin
1381  if (makeoldJEThadhistos) {
1382  ieta=GetEtaBin(deta); // bin of eta
1383  if(ieta<0) continue; // double check we don't have a negative array index
1384  iptjet=GetpTjetBin(jet->Pt()); // bin of jet pt
1385  if(iptjet<0) continue; // double check we don't have a negative array index
1386  }
1387 
1388  // dPHI between jet and hadron
1389  Double_t dphijh = RelativePhi(jet->Phi(), track->Phi()); // angle between jet and hadron
1390 
1391  // fill some jet-hadron histo's
1392  if (makeoldJEThadhistos) fHistJetH[centbin][iptjet][ieta]->Fill(dphijh,track->Pt()); // fill jet-hadron dPHI--track pT distribution
1393  if(makeQAhistos) fHistJetHEtaPhi->Fill(deta,dphijh); // fill jet-hadron eta--phi distribution
1394  fHistJetHaddPHI->Fill(dphijh);
1395  if(passedTTcut){
1396  if (makeoldJEThadhistos) fHistJetHTT[centbin][iptjet][ieta]->Fill(dphijh,track->Pt());
1397  }
1398 
1399  // does our max track or cluster pass the bias?
1400  if ((jet->MaxTrackPt()>fTrkBias) || (jet->MaxClusterPt()>fClusBias)){
1401  // set up and fill Jet-Hadron THnSparse
1402  if(fcorrJetPt){
1403  Double_t triggerEntries[9] = {fCent, jetPtLocal, track->Pt(), deta, dphijh, dEP, zVtx, trCharge, leadjet};
1404  if(fDoEventMixing) {
1405  if(fReduceStatsCent > 0) {
1406  if(cbin == fReduceStatsCent) fhnJH->Fill(triggerEntries, 1.0/trefficiency); // fill Sparse Histo with trigger entries
1407  } else fhnJH->Fill(triggerEntries, 1.0/trefficiency); // fill Sparse Histo with trigger entries
1408  }
1409  } else { // don't correct jet pt
1410  Double_t triggerEntries[9] = {fCent, jet->Pt(), track->Pt(), deta, dphijh, dEP, zVtx, trCharge, leadjet};
1411  if(fDoEventMixing) {
1412  if(fReduceStatsCent > 0) {
1413  if(cbin == fReduceStatsCent) fhnJH->Fill(triggerEntries, 1.0/trefficiency); // fill Sparse Histo with trigger entries
1414  } else fhnJH->Fill(triggerEntries, 1.0/trefficiency); // fill Sparse Histo with trigger entries
1415  }
1416  }
1417 
1418  // fill histo's
1419  if(makeQAhistos) fHistSEphieta->Fill(dphijh, deta); // single event distribution
1420  if(makeoldJEThadhistos) fHistJetHBias[centbin][iptjet][ieta]->Fill(dphijh,track->Pt());
1421 
1422  if (makeBIAShistos) {
1423  fHistJetHaddPhiBias->Fill(dphijh);
1424 
1425  // in plane and out of plane histo's
1426  if( dEP>0 && dEP<=(TMath::Pi()/6) ){
1427  // we are IN plane
1428  fHistJetHaddPhiINBias->Fill(dphijh);
1429  }else if( dEP>(TMath::Pi()/3) && dEP<=(TMath::Pi()/2) ){
1430  // we are OUT of PLANE
1431  fHistJetHaddPhiOUTBias->Fill(dphijh);
1432  }else if( dEP>(TMath::Pi()/6) && dEP<=(TMath::Pi()/3) ){
1433  // we are in middle of plane
1434  fHistJetHaddPhiMIDBias->Fill(dphijh);
1435  }
1436  } // BIAS histos switch
1437 
1438  } // end of check maxtrackpt>ftrackbias or maxclusterpt>fclustbias
1439 
1440  // **************************************************************************************************************
1441  // *********************************** PID **********************************************************************
1442  // **************************************************************************************************************
1443  if(doPIDtrackBIAS){
1444  //if(ptmax < fTrkBias) continue; // force PID to happen when max track pt > 5.0 GeV
1445  if(leadhadronPT < fTrkBias) continue; // force PID to happen when lead hadron pt > 5.0 GeV
1446  }
1447 
1448  // some variables for PID
1449  Double_t pt = -999, dEdx = -999, ITSsig = -999, TOFsig = -999, charge = -999;
1450 
1451  // nSigma of particles in TPC, TOF, and ITS
1452  Double_t nSigmaPion_TPC, nSigmaProton_TPC, nSigmaKaon_TPC;
1453  Double_t nSigmaPion_TOF, nSigmaProton_TOF, nSigmaKaon_TOF;
1454  Double_t nSigmaPion_ITS, nSigmaProton_ITS, nSigmaKaon_ITS;
1455 
1456  if(fPIDResponse) fHistTPCdEdX->Fill(track->Pt(), track->GetTPCsignal());
1457 
1458  if(doPID && (fPIDResponse) && ((jet->MaxTrackPt()>fTrkBias) || (jet->MaxClusterPt()>fClusBias)) ){
1459  // get parameters of track
1460  charge = track->Charge(); // charge of track
1461  pt = track->Pt(); // pT of track
1462 
1463  //if (!fPIDResponse) return kTRUE; // just return, maybe put at beginning
1464 
1465  fHistEventQA->Fill(13); // check for AliVEvent and fPIDresponse objects
1466 
1467  // get detector signals
1468  dEdx = track->GetTPCsignal();
1469  ITSsig = track->GetITSsignal();
1470  TOFsig = track->GetTOFsignal();
1471 
1472  // TPC nSigma's
1473  nSigmaPion_TPC = fPIDResponse->NumberOfSigmasTPC(track,AliPID::kPion);
1474  nSigmaKaon_TPC = fPIDResponse->NumberOfSigmasTPC(track,AliPID::kKaon);
1475  nSigmaProton_TPC = fPIDResponse->NumberOfSigmasTPC(track,AliPID::kProton);
1476 
1477  // TOF nSigma's
1478  nSigmaPion_TOF = fPIDResponse->NumberOfSigmasTOF(track,AliPID::kPion);
1479  nSigmaKaon_TOF = fPIDResponse->NumberOfSigmasTOF(track,AliPID::kKaon);
1480  nSigmaProton_TOF = fPIDResponse->NumberOfSigmasTOF(track,AliPID::kProton);
1481 
1482  // ITS nSigma's
1483  nSigmaPion_ITS = fPIDResponse->NumberOfSigmasITS(track,AliPID::kPion);
1484  nSigmaKaon_ITS = fPIDResponse->NumberOfSigmasITS(track,AliPID::kKaon);
1485  nSigmaProton_ITS = fPIDResponse->NumberOfSigmasITS(track,AliPID::kProton);
1486 
1487  // fill detector signal histograms
1488  //if (makeQAhistos) fHistTPCdEdX->Fill(pt, dEdx); // temp out
1489  if (makeQAhistos) fHistITSsignal->Fill(pt, ITSsig);
1490  //if (makeQAhistos) fHistTOFsignal->Fill(pt, TOFsig);
1491 
1492  // Tests to PID pions, kaons, and protons, (default is undentified tracks)
1493  //Double_t nPIDtpc = 0, nPIDits = 0, nPIDtof = 0;
1494  Double_t nPID = -99;
1495 
1496  // check track has pT < 0.900 GeV - use TPC pid
1497  if (pt<0.900 && dEdx>0) {
1498  //nPIDtpc = 4;
1499  nPID = 1;
1500 
1501  // PION check - TPC
1502  if (TMath::Abs(nSigmaPion_TPC)<2 && TMath::Abs(nSigmaKaon_TPC)>2 && TMath::Abs(nSigmaProton_TPC)>2 ){
1503  isPItpc = kTRUE;
1504  //nPIDtpc = 1;
1505  nPID=2;
1506  }else isPItpc = kFALSE;
1507 
1508  // KAON check - TPC
1509  if (TMath::Abs(nSigmaKaon_TPC)<2 && TMath::Abs(nSigmaPion_TPC)>3 && TMath::Abs(nSigmaProton_TPC)>2 ){
1510  isKtpc = kTRUE;
1511  //nPIDtpc = 2;
1512  nPID=3;
1513  }else isKtpc = kFALSE;
1514 
1515  // PROTON check - TPC
1516  if (TMath::Abs(nSigmaProton_TPC)<2 && TMath::Abs(nSigmaPion_TPC)>3 && TMath::Abs(nSigmaKaon_TPC)>2 ){
1517  isPtpc = kTRUE;
1518  //nPIDtpc = 3;
1519  nPID=4;
1520  }else isPtpc = kFALSE;
1521  } // cut on track pT for TPC
1522 
1523  // check track has pT < 0.500 GeV - use ITS pid
1524  if (pt<0.500 && ITSsig>0) {
1525  //nPIDits = 4;
1526  nPID = 5;
1527 
1528  // PION check - ITS
1529  if (TMath::Abs(nSigmaPion_ITS)<2 && TMath::Abs(nSigmaKaon_ITS)>2 && TMath::Abs(nSigmaProton_ITS)>2 ){
1530  isPIits = kTRUE;
1531  //nPIDits = 1;
1532  nPID=6;
1533  }else isPIits = kFALSE;
1534 
1535  // KAON check - ITS
1536  if (TMath::Abs(nSigmaKaon_ITS)<2 && TMath::Abs(nSigmaPion_ITS)>3 && TMath::Abs(nSigmaProton_ITS)>2 ){
1537  isKits = kTRUE;
1538  //nPIDits = 2;
1539  nPID=7;
1540  }else isKits = kFALSE;
1541 
1542  // PROTON check - ITS
1543  if (TMath::Abs(nSigmaProton_ITS)<2 && TMath::Abs(nSigmaPion_ITS)>3 && TMath::Abs(nSigmaKaon_ITS)>2 ){
1544  isPits = kTRUE;
1545  //nPIDits = 3;
1546  nPID=8;
1547  }else isPits = kFALSE;
1548  } // cut on track pT for ITS
1549 
1550  // check track has 0.900 GeV < pT < 2.500 GeV - use TOF pid
1551  if (pt>0.900 && pt<2.500 && TOFsig>0) {
1552  //nPIDtof = 4;
1553  nPID = 9;
1554 
1555  // PION check - TOF
1556  if (TMath::Abs(nSigmaPion_TOF)<2 && TMath::Abs(nSigmaKaon_TOF)>2 && TMath::Abs(nSigmaProton_TOF)>2 ){
1557  isPItof = kTRUE;
1558  //nPIDtof = 1;
1559  nPID=10;
1560  }else isPItof = kFALSE;
1561 
1562  // KAON check - TOF
1563  if (TMath::Abs(nSigmaKaon_TOF)<2 && TMath::Abs(nSigmaPion_TOF)>3 && TMath::Abs(nSigmaProton_TOF)>2 ){
1564  isKtof = kTRUE;
1565  //nPIDtof = 2;
1566  nPID=11;
1567  }else isKtof = kFALSE;
1568 
1569  // PROTON check - TOF
1570  if (TMath::Abs(nSigmaProton_TOF)<2 && TMath::Abs(nSigmaPion_TOF)>3 && TMath::Abs(nSigmaKaon_TOF)>2 ){
1571  isPtof = kTRUE;
1572  //nPIDtof = 3;
1573  nPID=12;
1574  }else isPtof = kFALSE;
1575  } // cut on track pT for TOF
1576 
1577  // if not tagged - unidentifed particle
1578  if (nPID == -99) nPID = 14;
1579 
1580  // fill PID tagged histo
1581  fHistPID->Fill(nPID);
1582 
1583  // TOF PID cuts: (needs testing)
1584  // Pion ID
1585  if((pt > 0.4) && (pt <=1.2) && (nSigmaPion_TOF >= -5.0) && (nSigmaPion_TOF <= 5.0)) nPID = 21;
1586  if((pt > 1.2) && (pt <=1.6) && (nSigmaPion_TOF >= -4.0) && (nSigmaPion_TOF <= 4.0)) nPID = 21;
1587  if((pt > 1.6) && (pt <=2.0) && (nSigmaPion_TOF >= -4.0) && (nSigmaPion_TOF <= 2.0)) nPID = 21;
1588  if((pt > 2.0) && (pt <=3.6) && (nSigmaPion_TOF >= -4.0) && (nSigmaPion_TOF <= 0.0)) nPID = 21;
1589  // Proton ID
1590  if((pt > 0.4) && (pt <=0.9) && (nSigmaProton_TOF >= -3.0) && (nSigmaPion_TOF <= 3.0)) nPID = 22;
1591  if((pt > 0.9) && (pt <=1.4) && (nSigmaProton_TOF >= -4.0) && (nSigmaPion_TOF <= 4.0)) nPID = 22;
1592  if((pt > 1.4) && (pt <=2.2) && (nSigmaProton_TOF >= -5.0) && (nSigmaPion_TOF <= 5.0)) nPID = 22;
1593  if((pt > 2.2) && (pt <=2.4) && (nSigmaProton_TOF >= -4.0) && (nSigmaPion_TOF <= 5.0)) nPID = 22;
1594  if((pt > 2.4) && (pt <=3.0) && (nSigmaProton_TOF >= -2.0) && (nSigmaPion_TOF <= 5.0)) nPID = 22;
1595  if((pt > 3.0) && (pt <=3.6) && (nSigmaProton_TOF >= 0.0) && (nSigmaPion_TOF <= 5.0)) nPID = 22;
1596  // Kaon ID
1597  if((pt > 0.4) && (pt <=0.8) && (nSigmaKaon_TOF >= -5.0) && (nSigmaPion_TOF <= 5.0)) nPID = 23;
1598  if((pt > 0.8) && (pt <=1.2) && (nSigmaKaon_TOF >= -4.0) && (nSigmaPion_TOF <= 4.0)) nPID = 23;
1599  if((pt > 1.2) && (pt <=1.7) && (nSigmaKaon_TOF >= -3.0) && (nSigmaPion_TOF <= 4.0)) nPID = 23;
1600  if((pt > 1.7) && (pt <=2.3) && (nSigmaKaon_TOF >= 0.0) && (nSigmaPion_TOF <= 5.0)) nPID = 23;
1601  if((pt > 2.3) && (pt <=2.6) && (nSigmaKaon_TOF >= 0.0) && (nSigmaPion_TOF <= 3.0)) nPID = 23;
1602  if((pt > 2.6) && (pt <=3.6) && (nSigmaKaon_TOF >= 0.0) && (nSigmaPion_TOF <= 2.0)) nPID = 23;
1603 
1604 
1605  // PID sparse getting filled
1606  if(fcorrJetPt){ // subtract background density
1607  if (allpidAXIS) { // FILL ALL axis
1608  Double_t pid_EntriesALL[19] = {fCent,pt,charge,deta,dphijh,leadjet,zVtx,dEP,jetPtLocal,
1609  nSigmaPion_TPC, nSigmaPion_TOF, // pion nSig values in TPC/TOF
1610  nPID, //nPIDtpc, nPIDits, nPIDtof, // PID label for each detector
1611  nSigmaProton_TPC, nSigmaKaon_TPC, // nSig in TPC
1612  nSigmaPion_ITS, nSigmaProton_ITS, nSigmaKaon_ITS, // nSig in ITS
1613  nSigmaProton_TOF, nSigmaKaon_TOF, // nSig in TOF
1614  }; //array for PID sparse
1615  if(fReduceStatsCent > 0) {
1616  if(cbin == fReduceStatsCent) fhnPID->Fill(pid_EntriesALL, 1.0/trefficiency);
1617  } else fhnPID->Fill(pid_EntriesALL, 1.0/trefficiency);
1618  } else {
1619  // PID sparse getting filled
1620  Double_t pid_Entries[12] = {fCent,pt,charge,deta,dphijh,leadjet,zVtx,dEP,jetPtLocal,
1621  nSigmaPion_TPC, nSigmaPion_TOF, // pion nSig values in TPC/TOF
1622  nPID //nPIDtpc, nPIDits, nPIDtof // PID label for each detector
1623  }; //array for PID sparse
1624  if(fReduceStatsCent > 0) {
1625  if(cbin == fReduceStatsCent) fhnPID->Fill(pid_Entries, 1.0/trefficiency);
1626  } else fhnPID->Fill(pid_Entries, 1.0/trefficiency); // fill Sparse histo of PID tracks
1627  } // minimal pid sparse filling
1628  } else { // don't correct jet pt by background density
1629  if (allpidAXIS) { // FILL ALL axis
1630  Double_t pid_EntriesALL[19] = {fCent,pt,charge,deta,dphijh,leadjet,zVtx,dEP,jetPt,
1631  nSigmaPion_TPC, nSigmaPion_TOF, // pion nSig values in TPC/TOF
1632  nPID, //nPIDtpc, nPIDits, nPIDtof, // PID label for each detector
1633  nSigmaProton_TPC, nSigmaKaon_TPC, // nSig in TPC
1634  nSigmaPion_ITS, nSigmaProton_ITS, nSigmaKaon_ITS, // nSig in ITS
1635  nSigmaProton_TOF, nSigmaKaon_TOF, // nSig in TOF
1636  }; //array for PID sparse
1637  if(fReduceStatsCent > 0) {
1638  if(cbin == fReduceStatsCent) fhnPID->Fill(pid_EntriesALL, 1.0/trefficiency);
1639  } else fhnPID->Fill(pid_EntriesALL, 1.0/trefficiency);
1640  } else {
1641  // PID sparse getting filled
1642  Double_t pid_Entries[12] = {fCent,pt,charge,deta,dphijh,leadjet,zVtx,dEP,jetPt,
1643  nSigmaPion_TPC, nSigmaPion_TOF, // pion nSig values in TPC/TOF
1644  nPID //nPIDtpc, nPIDits, nPIDtof // PID label for each detector
1645  }; //array for PID sparse
1646  if(fReduceStatsCent > 0) {
1647  if(cbin == fReduceStatsCent) fhnPID->Fill(pid_Entries, 1.0/trefficiency);
1648  } else fhnPID->Fill(pid_Entries, 1.0/trefficiency);
1649  } // minimal pid sparse filling
1650  } // don't correct jet pet
1651  } // end of doPID check
1652 
1653  // get track pt bin
1654  Int_t itrackpt = -500; // initialize track pT bin
1655  itrackpt = GetpTtrackBin(track->Pt());
1656 
1657  // all tracks: jet hadron correlations in hadron pt bins
1658  if(makeextraCORRhistos) fHistJetHadbindPhi[itrackpt]->Fill(dphijh);
1659 
1660  // in plane and out of plane jet-hadron histo's
1661  if( dEP>0 && dEP<=(TMath::Pi()/6) ){
1662  // we are IN plane
1663  if(makeextraCORRhistos) fHistJetHaddPhiINcent[centbin]->Fill(dphijh);
1664  fHistJetHaddPhiIN->Fill(dphijh);
1665  if(makeextraCORRhistos) fHistJetHadbindPhiIN[itrackpt]->Fill(dphijh);
1666  }else if( dEP>(TMath::Pi()/3) && dEP<=(TMath::Pi()/2) ){
1667  // we are OUT of PLANE
1668  if(makeextraCORRhistos) fHistJetHaddPhiOUTcent[centbin]->Fill(dphijh);
1669  fHistJetHaddPhiOUT->Fill(dphijh);
1670  if(makeextraCORRhistos) fHistJetHadbindPhiOUT[itrackpt]->Fill(dphijh);
1671  }else if( dEP>(TMath::Pi()/6) && dEP<=(TMath::Pi()/3) ){
1672  // we are in the middle of plane
1673  if(makeextraCORRhistos) fHistJetHaddPhiMIDcent[centbin]->Fill(dphijh);
1674  fHistJetHaddPhiMID->Fill(dphijh);
1675  if(makeextraCORRhistos) fHistJetHadbindPhiMID[itrackpt]->Fill(dphijh);
1676  }
1677  } // loop over tracks found in event with highest JET pT > 10.0 GeV (change)
1678  } // jet pt cut
1679  } // jet loop
1680 
1681  fHistEventQA->Fill(14); // events right before event mixing
1682 
1683 // ***************************************************************************************************************
1684 // ******************************** Event MIXING *****************************************************************
1685 // ***************************************************************************************************************
1686 
1687  // initialize object array of cloned picotracks
1688  TObjArray* tracksClone = 0x0;
1689 
1690  // PbPb collisions - create cloned picotracks
1691  //if(GetBeamType() == 1) tracksClone = CloneAndReduceTrackList(tracks); // TEST
1692 
1693  //Prepare to do event mixing
1694  if(fDoEventMixing>0){
1695  // event mixing
1696 
1697  // 1. First get an event pool corresponding in mult (cent) and
1698  // zvertex to the current event. Once initialized, the pool
1699  // should contain nMix (reduced) events. This routine does not
1700  // pre-scan the chain. The first several events of every chain
1701  // will be skipped until the needed pools are filled to the
1702  // specified depth. If the pool categories are not too rare, this
1703  // should not be a problem. If they are rare, you could lose
1704  // statistics.
1705 
1706  // 2. Collect the whole pool's content of tracks into one TObjArray
1707  // (bgTracks), which is effectively a single background super-event.
1708 
1709  // 3. The reduced and bgTracks arrays must both be passed into
1710  // FillCorrelations(). Also nMix should be passed in, so a weight
1711  // of 1./nMix can be applied.
1712 
1713  // mix jets from triggered events with tracks from MB events
1714  // get the trigger bit, need to change trigger bits between different runs
1715  UInt_t trigger = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
1716  // if event was not selected (triggered) for any reseason (should never happen) then return
1717  if (trigger==0) return kTRUE;
1718 
1719  // initialize event pools
1720  AliEventPool* pool = 0x0;
1721  AliEventPool* poolpp = 0x0;
1722  Double_t Ntrks = -999;
1723 
1724  // pp collisions - get event pool
1725  if(GetBeamType() == 0) {
1726  Ntrks=(Double_t)Ntracks*1.0;
1727  //cout<<"Test.. Ntrks: "<<fPoolMgr->GetEventPool(Ntrks);
1728  poolpp = fPoolMgr->GetEventPool(Ntrks, zVtx); // for pp
1729  }
1730 
1731  // PbPb collisions - get event pool
1732  if(GetBeamType() == 1) pool = fPoolMgr->GetEventPool(fCent, zVtx); // for PbPb? fcent
1733 
1734  // if we don't have a pool, return
1735  if (!pool && !poolpp){
1736  if(GetBeamType() == 1) AliFatal(Form("No pool found for centrality = %f, zVtx = %f", fCent, zVtx));
1737  if(GetBeamType() == 0) AliFatal(Form("No pool found for multiplicity = %f, zVtx = %f", Ntrks, zVtx));
1738  return kTRUE;
1739  }
1740 
1741  fHistEventQA->Fill(15); // mixed events cases that have pool
1742 
1743  // initialize background tracks array
1744  TObjArray* bgTracks;
1745 
1746  // next line might not apply for PbPb collisions
1747  // use only jets from EMCal-triggered events (for lhc11a use AliVEvent::kEMC1)
1748  //check for a trigger jet
1749  // fmixingtrack/10 ??
1750  if(GetBeamType() == 1) if(trigger && fTriggerEventType) { //kEMCEJE)) {
1751  if (pool->IsReady() || pool->NTracksInPool() > fNMIXtracks || pool->GetCurrentNEvents() >= fNMIXevents) {
1752 
1753  // loop over jets (passing cuts?)
1754  for (Int_t ijet = 0; ijet < Njets; ijet++) {
1755  Double_t leadjet=0;
1756  if (ijet==ijethi) leadjet=1;
1757 
1758  // get jet object
1759  AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(ijet));
1760  if (!jet) continue;
1761 
1762  // (should probably be higher..., but makes a cut on jet pT)
1763  if (jet->Pt()<0.1) continue;
1764  if (!AcceptMyJet(jet)) continue;
1765 
1766  Double_t jetPtLocalmix = -500; // initialize corr jet pt
1767  if(GetBeamType() == 1) {
1768  fLocalRhoVal = fLocalRho->GetLocalVal(jet->Phi(), fJetRad); //GetJetRadius(0)); // get local rho value
1769  jetPtLocalmix = jet->Pt()-jet->Area()*fLocalRhoVal; // corrected pT of jet using Rho modulated for V2 and V3
1770  }
1771 
1772  fHistEventQA->Fill(16); // event mixing jets
1773 
1774  // set cut to do event mixing only if we have a jet meeting our pt threshold (bias applied below)
1775  if (jet->Pt()<fJetPtcut) continue;
1776 
1777  // get number of current events in pool
1778  Int_t nMix = pool->GetCurrentNEvents(); // how many particles in pool to mix
1779 
1780  // Fill for biased jet triggers only
1781  if ((jet->MaxTrackPt()>fTrkBias) || (jet->MaxClusterPt()>fClusBias)) { // && jet->Pt() > fJetPtcut) {
1782  // Fill mixed-event histos here
1783  for (Int_t jMix=0; jMix<nMix; jMix++) {
1784  fHistEventQA->Fill(17); // event mixing nMix
1785 
1786  // get jMix'th event
1787  bgTracks = pool->GetEvent(jMix);
1788  const Int_t Nbgtrks = bgTracks->GetEntries();
1789  for(Int_t ibg=0; ibg<Nbgtrks; ibg++) {
1790  AliPicoTrack *part = static_cast<AliPicoTrack*>(bgTracks->At(ibg));
1791  if(!part) continue;
1792  if(TMath::Abs(part->Eta())>0.9) continue;
1793  if(part->Pt()<0.15) continue;
1794 
1795  Double_t DEta = part->Eta()-jet->Eta(); // difference in eta
1796  Double_t DPhi = RelativePhi(jet->Phi(),part->Phi()); // difference in phi
1797  Double_t dEP = RelativeEPJET(jet->Phi(),fEPV0); // difference between jet and EP
1798  Double_t mixcharge = part->Charge();
1799  //Double_t DR=TMath::Sqrt(DPhi*DPhi+DEta*DEta); // difference in R
1800 
1801  // calculate single particle tracking efficiency of mixed events for correlations
1802  Double_t mixefficiency = -999;
1803  mixefficiency = EffCorrection(part->Eta(), part->Pt(), fDoEffCorr);
1804 
1805  // create / fill mixed event sparse
1806  if(fcorrJetPt) {
1807  Double_t triggerEntries[9] = {fCent,jetPtLocalmix,part->Pt(),DEta,DPhi,dEP,zVtx, mixcharge, leadjet}; //array for ME sparse
1808  if(fReduceStatsCent > 0) {
1809  if(cbin == fReduceStatsCent) fhnMixedEvents->Fill(triggerEntries,1./(nMix*mixefficiency)); // fill Sparse histo of mixed events
1810  } else fhnMixedEvents->Fill(triggerEntries,1./(nMix*mixefficiency)); // fill Sparse histo of mixed events
1811  } else { // don't correct jet pt
1812  Double_t triggerEntries[9] = {fCent,jet->Pt(),part->Pt(),DEta,DPhi,dEP,zVtx, mixcharge, leadjet}; //array for ME sparse
1813  if(fReduceStatsCent > 0) {
1814  if(cbin == fReduceStatsCent) fhnMixedEvents->Fill(triggerEntries,1./(nMix*mixefficiency)); // fill Sparse histo of mixed events
1815  } else fhnMixedEvents->Fill(triggerEntries,1./(nMix*mixefficiency)); // fill Sparse histo of mixed events
1816  }
1817 
1818  fHistEventQA->Fill(18); // event mixing - nbgtracks
1819  if(makeQAhistos) fHistMEphieta->Fill(DPhi,DEta, 1./(nMix*mixefficiency));
1820  } // end of background track loop
1821  } // end of filling mixed-event histo's
1822  } // end of check for biased jet triggers
1823  } // end of jet loop
1824  } // end of check for pool being ready
1825  } // end EMC triggered loop
1826 
1827 //=============================================================================================================
1828 
1829  // use only jets from EMCal-triggered events (for lhc11a use AliVEvent::kEMC1)
1831  // pp collisions
1832  if(GetBeamType() == 0) if(trigger && fTriggerEventType) { //kEMC1)) {
1833  if (poolpp->IsReady() || poolpp->NTracksInPool() > fNMIXtracks || poolpp->GetCurrentNEvents() >= fNMIXevents) {
1834 
1835  // loop over jets (passing cuts?)
1836  for (Int_t ijet = 0; ijet < Njets; ijet++) {
1837  Double_t leadjet=0;
1838  if (ijet==ijethi) leadjet=1;
1839 
1840  // get jet object
1841  AliEmcalJet *jet = static_cast<AliEmcalJet*>(jets->At(ijet));
1842  if (!jet) continue;
1843 
1844  // (should probably be higher..., but makes a cut on jet pT)
1845  if (jet->Pt()<0.1) continue;
1846  if (!AcceptMyJet(jet)) continue;
1847 
1848  fHistEventQA->Fill(16); // event mixing jets
1849 
1850  // set cut to do event mixing only if we have a jet meeting our pt threshold (bias applied below)
1851  if (jet->Pt()<fJetPtcut) continue;
1852 
1853  // get number of current events in pool
1854  Int_t nMix = poolpp->GetCurrentNEvents(); // how many particles in pool to mix
1855 
1856  // Fill for biased jet triggers only
1857  if ((jet->MaxTrackPt()>fTrkBias) || (jet->MaxClusterPt()>fClusBias)) { // && jet->Pt() > fJetPtcut) {
1858  // Fill mixed-event histos here
1859  for (Int_t jMix=0; jMix<nMix; jMix++) {
1860  fHistEventQA->Fill(17); // event mixing nMix
1861 
1862  // get jMix'th event
1863  bgTracks = poolpp->GetEvent(jMix);
1864  const Int_t Nbgtrks = bgTracks->GetEntries();
1865  for(Int_t ibg=0; ibg<Nbgtrks; ibg++) {
1866  AliPicoTrack *part = static_cast<AliPicoTrack*>(bgTracks->At(ibg));
1867  if(!part) continue;
1868  if(TMath::Abs(part->Eta())>0.9) continue;
1869  if(part->Pt()<0.15) continue;
1870 
1871  Double_t DEta = part->Eta()-jet->Eta(); // difference in eta
1872  Double_t DPhi = RelativePhi(jet->Phi(),part->Phi()); // difference in phi
1873  Double_t dEP = RelativeEPJET(jet->Phi(),fEPV0); // difference between jet and EP
1874  Double_t mixcharge = part->Charge();
1875  //Double_t DR=TMath::Sqrt(DPhi*DPhi+DEta*DEta); // difference in R
1876 
1877  // create / fill mixed event sparse
1878  Double_t triggerEntries[9] = {fCent,jet->Pt(),part->Pt(),DEta,DPhi,dEP,zVtx, mixcharge, leadjet}; //array for ME sparse
1879  fhnMixedEvents->Fill(triggerEntries,1./nMix); // fill Sparse histo of mixed events
1880 
1881  fHistEventQA->Fill(18); // event mixing - nbgtracks
1882  if(makeextraCORRhistos) fHistMEphieta->Fill(DPhi,DEta, 1./nMix);
1883  } // end of background track loop
1884  } // end of filling mixed-event histo's
1885  } // end of check for biased jet triggers
1886  } // end of jet loop
1887  } // end of check for pool being ready
1888  } //end EMC triggered loop
1889 
1890  // pp collisions
1891  if(GetBeamType() == 0) {
1892 
1893  // use only tracks from MB events (for lhc11a use AliVEvent::kMB) //kAnyINT as test
1894  if(trigger && fMixingEventType) { //kMB) {
1895 
1896  // create a list of reduced objects. This speeds up processing and reduces memory consumption for the event pool
1897  tracksClone = CloneAndReduceTrackList(tracks);
1898 
1899  // update pool if jet in event or not
1900  poolpp->UpdatePool(tracksClone);
1901 
1902  } // check on track from MB events
1903  }
1904 
1905  // PbPb collisions
1906  if(GetBeamType() == 1) {
1907 
1908  // use only tracks from MB and Central and Semi-Central events
1909  if(trigger && fMixingEventType) { //kMB) {
1910 
1911  // create a list of reduced objects. This speeds up processing and reduces memory consumption for the event pool
1912  tracksClone = CloneAndReduceTrackList(tracks);
1913 
1914  // update pool if jet in event or not
1915  pool->UpdatePool(tracksClone);
1916 
1917  } // MB and Central and Semi-Central events
1918  } // PbPb collisions
1919  } // end of event mixing
1920 
1921  // print some stats on the event
1922  event++;
1923  fHistEventQA->Fill(19); // events making it to end
1924 
1925  // fill Event Trigger QA (after cuts)
1927 
1928  if (doComments) {
1929  cout<<"Event #: "<<event<<" Jet Radius: "<<fJetRad<<" Constituent Pt Cut: "<<fConstituentCut<<endl;
1930  cout<<"# of jets: "<<Njets<<" NjetAcc: "<<NjetAcc<<" Highest jet pt: "<<highestjetpt<<" leading hadron pt: "<<leadhadronPT<<endl;
1931  cout<<"# tracks: "<<Ntracks<<" NtrackAcc: "<<NtrackAcc<<" Highest track pt: "<<ptmax<<endl;
1932  cout<<"maxcluster pt = "<<maxclusterpt<<" maxtrack pt = "<<maxtrackpt<<endl;
1933  cout<<" =============================================== "<<endl;
1934  }
1935 
1936  return kTRUE; // used when the function is of type bool
1937 } // end of RUN
1938 
1939 //________________________________________________________________________
1941 { // Get centrality bin.
1942  Int_t centbin = -1;
1943  if (cent>=0 && cent<10) centbin = 0;
1944  else if (cent>=10 && cent<20) centbin = 1;
1945  else if (cent>=20 && cent<30) centbin = 2;
1946  else if (cent>=30 && cent<40) centbin = 3;
1947  else if (cent>=40 && cent<50) centbin = 4;
1948  else if (cent>=50 && cent<90) centbin = 5;
1949 
1950  return centbin;
1951 }
1952 
1953 //________________________________________________________________________
1954 Double_t AliAnalysisTaskEmcalJetHadEPpid::RelativePhi(Double_t mphi,Double_t vphi) const
1955 { // function to calculate relative PHI
1956  double dphi = mphi-vphi;
1957 
1958  // set dphi to operate on adjusted scale
1959  if(dphi<-0.5*TMath::Pi()) dphi+=2.*TMath::Pi();
1960  if(dphi>3./2.*TMath::Pi()) dphi-=2.*TMath::Pi();
1961 
1962  // test
1963  if( dphi < -1.*TMath::Pi()/2 || dphi > 3.*TMath::Pi()/2 )
1964  AliWarning(Form("%s: dPHI not in range [-0.5*Pi, 1.5*Pi]!", GetName()));
1965 
1966  return dphi; // dphi in [-0.5Pi, 1.5Pi]
1967 }
1968 
1969 //_________________________________________________________________________
1970 Double_t AliAnalysisTaskEmcalJetHadEPpid::RelativeEPJET(Double_t jetAng, Double_t EPAng) const
1971 { // function to calculate angle between jet and EP in the 1st quadrant (0,Pi/2)
1972  Double_t dphi = (EPAng - jetAng);
1973 
1974  // ran into trouble with a few dEP<-Pi so trying this...
1975  if( dphi<-1*TMath::Pi() ){
1976  dphi = dphi + 1*TMath::Pi();
1977  } // this assumes we are doing full jets currently
1978 
1979  if( (dphi>0) && (dphi<1*TMath::Pi()/2) ){
1980  // Do nothing! we are in quadrant 1
1981  }else if( (dphi>1*TMath::Pi()/2) && (dphi<1*TMath::Pi()) ){
1982  dphi = 1*TMath::Pi() - dphi;
1983  }else if( (dphi<0) && (dphi>-1*TMath::Pi()/2) ){
1984  dphi = fabs(dphi);
1985  }else if( (dphi<-1*TMath::Pi()/2) && (dphi>-1*TMath::Pi()) ){
1986  dphi = dphi + 1*TMath::Pi();
1987  }
1988 
1989  // test
1990  if( dphi < 0 || dphi > TMath::Pi()/2 )
1991  AliWarning(Form("%s: dPHI not in range [0, 0.5*Pi]!", GetName()));
1992 
1993  return dphi; // dphi in [0, Pi/2]
1994 }
1995 
1996 //Int_t ieta=GetEtaBin(deta);
1997 //________________________________________________________________________
1999 {
2000  // Get eta bin for histos.
2001  Int_t etabin = -1;
2002  if (TMath::Abs(eta)<=0.4) etabin = 0;
2003  else if (TMath::Abs(eta)>0.4 && TMath::Abs(eta)<0.8) etabin = 1;
2004  else if (TMath::Abs(eta)>=0.8) etabin = 2;
2005 
2006  return etabin;
2007 } // end of get-eta-bin
2008 
2009 //________________________________________________________________________
2011 {
2012  // Get jet pt bin for histos.
2013  Int_t ptbin = -1;
2014  if (pt>=15 && pt<20) ptbin = 0;
2015  else if (pt>=20 && pt<25) ptbin = 1;
2016  else if (pt>=25 && pt<40) ptbin = 2;
2017  else if (pt>=40 && pt<60) ptbin = 3;
2018  else if (pt>=60) ptbin = 4;
2019 
2020  return ptbin;
2021 } // end of get-jet-pt-bin
2022 
2023 //________________________________________________________________________
2025 {
2026  // May need to update bins for future runs... (testing locally)
2027 
2028  // Get track pt bin for histos.
2029  Int_t ptbin = -1;
2030  if (pt < 0.5) ptbin = 0;
2031  else if (pt>=0.5 && pt<1.0) ptbin = 1;
2032  else if (pt>=1.0 && pt<1.5) ptbin = 2;
2033  else if (pt>=1.5 && pt<2.0) ptbin = 3;
2034  else if (pt>=2.0 && pt<2.5) ptbin = 4;
2035  else if (pt>=2.5 && pt<3.0) ptbin = 5;
2036  else if (pt>=3.0 && pt<4.0) ptbin = 6;
2037  else if (pt>=4.0 && pt<5.0) ptbin = 7;
2038  else if (pt>=5.0) ptbin = 8;
2039 
2040  return ptbin;
2041 } // end of get-jet-pt-bin
2042 
2043 //___________________________________________________________________________
2045 {
2046  // get z-vertex bin for histo.
2047  int zVbin= -1;
2048  if (zVtx>=-10 && zVtx<-8) zVbin = 0;
2049  else if (zVtx>=-8 && zVtx<-6) zVbin = 1;
2050  else if (zVtx>=-6 && zVtx<-4) zVbin = 2;
2051  else if (zVtx>=-4 && zVtx<-2) zVbin = 3;
2052  else if (zVtx>=-2 && zVtx<0) zVbin = 4;
2053  else if (zVtx>=0 && zVtx<2) zVbin = 5;
2054  else if (zVtx>=2 && zVtx<4) zVbin = 6;
2055  else if (zVtx>=4 && zVtx<6) zVbin = 7;
2056  else if (zVtx>=6 && zVtx<8) zVbin = 8;
2057  else if (zVtx>=8 && zVtx<10) zVbin = 9;
2058  else zVbin = 10;
2059 
2060  return zVbin;
2061 } // end of get z-vertex bin
2062 
2063 //______________________________________________________________________
2064 THnSparse* AliAnalysisTaskEmcalJetHadEPpid::NewTHnSparseF(const char* name, UInt_t entries)
2065 {
2066  // generate new THnSparseF, axes are defined in GetDimParams()
2067  Int_t count = 0;
2068  UInt_t tmp = entries;
2069  while(tmp!=0){
2070  count++;
2071  tmp = tmp &~ -tmp; // clear lowest bit
2072  }
2073 
2074  TString hnTitle(name);
2075  const Int_t dim = count;
2076  Int_t nbins[dim];
2077  Double_t xmin[dim];
2078  Double_t xmax[dim];
2079 
2080  Int_t i=0;
2081  Int_t c=0;
2082  while(c<dim && i<32){
2083  if(entries&(1<<i)){
2084  TString label("");
2085  GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
2086  hnTitle += Form(";%s",label.Data());
2087  c++;
2088  }
2089 
2090  i++;
2091  }
2092  hnTitle += ";";
2093 
2094  return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
2095 } // end of NewTHnSparseF
2096 
2097 void AliAnalysisTaskEmcalJetHadEPpid::GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
2098 {
2099  // stores label and binning of axis for THnSparse
2100  const Double_t pi = TMath::Pi();
2101 
2102  switch(iEntry){
2103 
2104  case 0:
2105  label = "V0 centrality (%)";
2106  nbins = 10;
2107  xmin = 0.;
2108  xmax = 100.;
2109  break;
2110 
2111  case 1:
2112  if(fcorrJetPt) { // correct jet pt
2113  label = "Jet Corrected p_{T}";
2114  nbins = 50;
2115  xmin = -50.;
2116  xmax = 200.;
2117  } else { // don't correct jet pt
2118  label = "Jet p_{T}";
2119  if(doaltPIDbinning) { // minimize bins (for PID)
2120  nbins = 20;
2121  xmin = 0.;
2122  xmax = 100.;
2123  } else { // don't minimize bins
2124  nbins = 50;
2125  xmin = 0.;
2126  xmax = 250.;
2127  }
2128  }
2129  break;
2130 
2131  case 2:
2132  label = "Track p_{T}";
2133  if(doaltPIDbinning) { // minimize bins (for PID)
2134  nbins = 100;
2135  xmin = 0.;
2136  xmax = 10.;
2137  } else { // don't minimize bins
2138  nbins = 80;
2139  xmin = 0.;
2140  xmax = 20.;
2141  }
2142  break;
2143 
2144  case 3:
2145  label = "Relative Eta";
2146  nbins = 56; // 48
2147  xmin = -1.4;
2148  xmax = 1.4;
2149  break;
2150 
2151  case 4:
2152  label = "Relative Phi";
2153  nbins = 72;
2154  xmin = -0.5*pi;
2155  xmax = 1.5*pi;
2156  break;
2157 
2158  case 5:
2159  label = "Relative angle of Jet and Reaction Plane";
2160  nbins = 3; // (12) 72
2161  xmin = 0;
2162  xmax = 0.5*pi;
2163  break;
2164 
2165  case 6:
2166  label = "z-vertex";
2167  nbins = 10;
2168  xmin = -10;
2169  xmax = 10;
2170  break;
2171 
2172  case 7:
2173  label = "track charge";
2174  nbins = 3;
2175  xmin = -1.5;
2176  xmax = 1.5;
2177  break;
2178 
2179  case 8:
2180  label = "leading jet";
2181  if(doaltPIDbinning) nbins = 1;
2182  else nbins = 3;
2183  xmin = -0.5;
2184  xmax = 2.5;
2185  break;
2186 
2187  case 9: // need to update
2188  label = "leading track";
2189  if(doaltPIDbinning) nbins = 1;
2190  else nbins = 10;
2191  xmin = 0;
2192  xmax = 50;
2193  break;
2194 
2195  } // end of switch
2196 } // end of getting dim-params
2197 
2198 //_________________________________________________
2199 // From CF event mixing code PhiCorrelations
2201 {
2202  // clones a track list by using AliPicoTrack which uses much less memory (used for event mixing)
2203  TObjArray* tracksClone = new TObjArray;
2204  tracksClone->SetOwner(kTRUE);
2205 
2206 // ===============================
2207 
2208 // cout << "RM Hybrid track : " << i << " " << particle->Pt() << endl;
2209 
2210  //cout << "nEntries " << tracks->GetEntriesFast() <<endl;
2211  for (Int_t i=0; i<tracksME->GetEntriesFast(); i++) { // AOD/general case
2212  AliVParticle* particle = (AliVParticle*) tracksME->At(i); // AOD/general case
2213  if(TMath::Abs(particle->Eta())>fTrkEta) continue;
2214  if(particle->Pt()<0.15)continue;
2215 
2216 /*
2217 // DON'T USE
2218  Double_t trackpt=particle->Pt(); // track pT
2219 
2220  Int_t trklabel=-1;
2221  trklabel=particle->GetLabel();
2222  //cout << "TRACK_LABEL: " << particle->GetLabel()<<endl;
2223 
2224  Int_t hadbin=-1;
2225  if(trackpt<0.5) hadbin=0;
2226  else if(trackpt<1) hadbin=1;
2227  else if(trackpt<2) hadbin=2;
2228  else if(trackpt<3) hadbin=3;
2229  else if(trackpt<5) hadbin=4;
2230  else if(trackpt<8) hadbin=5;
2231  else if(trackpt<20) hadbin=6;
2232 // end of DON'T USE
2233 
2234 //feb10 comment out
2235  if(hadbin>-1 && trklabel>-1 && trklabel <3) fHistTrackEtaPhi[trklabel][hadbin]->Fill(particle->Eta(),particle->Phi());
2236  if(hadbin>-1) fHistTrackEtaPhi[3][hadbin]->Fill(particle->Eta(),particle->Phi());
2237 
2238  if(hadbin>-1) fHistTrackEtaPhi[hadbin]->Fill(particle->Eta(),particle->Phi()); // TEST
2239 */
2240 
2241  tracksClone->Add(new AliPicoTrack(particle->Pt(), particle->Eta(), particle->Phi(), particle->Charge(), 0, 0, 0, 0));
2242  } // end of looping through tracks
2243 
2244  return tracksClone;
2245 }
2246 
2247 //____________________________________________________________________________________________
2248 THnSparse* AliAnalysisTaskEmcalJetHadEPpid::NewTHnSparseFPID(const char* name, UInt_t entries)
2249 {
2250  // generate new THnSparseF PID, axes are defined in GetDimParams()
2251  Int_t count = 0;
2252  UInt_t tmp = entries;
2253  while(tmp!=0){
2254  count++;
2255  tmp = tmp &~ -tmp; // clear lowest bit
2256  }
2257 
2258  TString hnTitle(name);
2259  const Int_t dim = count;
2260  Int_t nbins[dim];
2261  Double_t xmin[dim];
2262  Double_t xmax[dim];
2263 
2264  Int_t i=0;
2265  Int_t c=0;
2266  while(c<dim && i<32){
2267  if(entries&(1<<i)){
2268  TString label("");
2269  GetDimParamsPID(i, label, nbins[c], xmin[c], xmax[c]);
2270  hnTitle += Form(";%s",label.Data());
2271  c++;
2272  }
2273 
2274  i++;
2275  }
2276  hnTitle += ";";
2277 
2278  return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
2279 } // end of NewTHnSparseF PID
2280 
2281 //________________________________________________________________________________
2282 void AliAnalysisTaskEmcalJetHadEPpid::GetDimParamsPID(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
2283 {
2284  // stores label and binning of axis for THnSparse
2285  const Double_t pi = TMath::Pi();
2286 
2287  switch(iEntry){
2288 
2289  case 0:
2290  label = "V0 centrality (%)";
2291  nbins = 10;
2292  xmin = 0.;
2293  xmax = 100.;
2294  break;
2295 
2296  case 1:
2297  label = "Track p_{T}";
2298  if(doaltPIDbinning) { // minimize bins
2299  nbins = 100;
2300  xmin = 0.;
2301  xmax = 10.;
2302  } else { // don't minmize
2303  nbins = 80;
2304  xmin = 0.;
2305  xmax = 20.;
2306  }
2307  break;
2308 
2309  case 2:
2310  label = "Charge of Track";
2311  nbins = 3;
2312  xmin = -1.5;
2313  xmax = 1.5;
2314  break;
2315 
2316  case 3:
2317  label = "Relative Eta of Track and Jet";
2318  nbins = 56; // 48
2319  xmin = -1.4;
2320  xmax = 1.4;
2321  break;
2322 
2323  case 4:
2324  label = "Relative Phi of Track and Jet";
2325  nbins = 72;
2326  xmin = -0.5*pi;
2327  xmax = 1.5*pi;
2328  break;
2329 
2330  case 5: // not used
2331  label = "leading jet";
2332  if(doaltPIDbinning) nbins = 1;
2333  else nbins = 3;
2334  xmin = -.5;
2335  xmax = 2.5;
2336  break;
2337 
2338  case 6:
2339  label = "Z-vertex";
2340  nbins = 10;
2341  xmin = -10.;
2342  xmax = 10.;
2343  break;
2344 
2345  case 7:
2346  label = "Relative angle: Jet and Reaction Plane";
2347  nbins = 3; // (12) 48
2348  xmin = 0.;
2349  xmax = 0.5*pi;
2350  break;
2351 
2352  case 8:
2353  if(fcorrJetPt) { // correct jet pt
2354  label = "Jet Corrected p_{T}";
2355  nbins = 50;
2356  xmin = -50.;
2357  xmax = 200.;
2358  } else { // don't correct jet pt
2359  label = "Jet p_{T}";
2360  if(doaltPIDbinning) { // minimize bins
2361  nbins = 20;
2362  xmin = 0.;
2363  xmax = 100.;
2364  } else { // don't minimize bins
2365  nbins = 50;
2366  xmin = 0.;
2367  xmax = 250.;
2368  }
2369  }
2370  break;
2371 
2372  case 9: // this may be temp, only using TOF
2373  label = "N-Sigma of pions in TPC";
2374  if(doaltPIDbinning) nbins = 1;
2375  else nbins = 200;
2376  xmin = -10.0;
2377  xmax = 10.0;
2378  break;
2379 
2380  case 10:
2381  label = "N-Sigma of pions in TOF";
2382  if(doaltPIDbinning) { // this may be temp, only using TOF
2383  nbins = 10;
2384  xmin = -5.;
2385  xmax = 5.;
2386  } else {
2387  nbins = 200;
2388  xmin = -10.;
2389  xmax = 10.;
2390  }
2391  break;
2392 
2393  case 11:
2394  label = "PID determination TPC 1-15, TOF 21-25";
2395  nbins = 25; //5;
2396  xmin = 0.5; //0.;
2397  xmax = 25.5; //5.;
2398  break;
2399 
2400 /*
2401  case 12:
2402  label = "ITS PID determination";
2403  nbins = 5;
2404  xmin = 0.;
2405  xmax = 5.;
2406  break;
2407 
2408  case 13:
2409  label = "TOF PID determination";
2410  nbins = 5;
2411  xmin = 0.;
2412  xmax = 5.;
2413  break;
2414 */
2415 
2416  case 12: // this may be temp, only using TOF
2417  label = "N-Sigma of protons in TPC";
2418  if(doaltPIDbinning) nbins = 1;
2419  else nbins = 200;
2420  xmin = -10.;
2421  xmax = 10.;
2422  break;
2423 
2424  case 13: // this may be temp, only using TOF
2425  label = "N-Sigma of kaons in TPC";
2426  if(doaltPIDbinning) nbins = 1;
2427  else nbins = 200;
2428  xmin = -10.;
2429  xmax = 10.;
2430  break;
2431 
2432  case 14: // this may be temp, only using TOF
2433  label = "N-Sigma of pions in ITS";
2434  if(doaltPIDbinning) nbins = 1;
2435  else nbins = 200;
2436  xmin = -10.0;
2437  xmax = 10.0;
2438  break;
2439 
2440  case 15: // this may be temp, only using TOF
2441  label = "N-Sigma of protons in ITS";
2442  if(doaltPIDbinning) nbins = 1;
2443  else nbins = 200;
2444  xmin = -10.;
2445  xmax = 10.;
2446  break;
2447 
2448  case 16: // this may be temp, only using TOF
2449  label = "N-Sigma of kaons in ITS";
2450  if(doaltPIDbinning) nbins = 1;
2451  else nbins = 200;
2452  xmin = -10.;
2453  xmax = 10.;
2454  break;
2455 
2456  case 17:
2457  label = "N-Sigma of protons in TOF";
2458  if(doaltPIDbinning) {
2459  nbins = 10;
2460  xmin = -5.;
2461  xmax = 5.;
2462  } else {
2463  nbins = 200;
2464  xmin = -10.;
2465  xmax = 10.;
2466  }
2467  break;
2468 
2469  case 18:
2470  label = "N-Sigma of kaons in TOF";
2471  if(doaltPIDbinning) {
2472  nbins = 10;
2473  xmin = -5.;
2474  xmax = 5.;
2475  } else {
2476  nbins = 200;
2477  xmin = -10.;
2478  xmax = 10.;
2479  }
2480  break;
2481 
2482  } // end of switch
2483 } // end of get dimension parameters PID
2484 
2486  cout<<"#########################"<<endl;
2487  cout<<"#### DONE RUNNING!!! ####"<<endl;
2488  cout<<"#########################"<<endl;
2489 } // end of terminate
2490 
2491 //________________________________________________________________________
2493  //applies all jet cuts except pt
2494  if ((jet->Phi()<fPhimin)||(jet->Phi()>fPhimax)) return 0;
2495  if ((jet->Eta()<fEtamin)||(jet->Eta()>fEtamax)) return 0;
2496  if (jet->Area()<fAreacut) return 0;
2497  // prevents 0 area jets from sneaking by when area cut == 0
2498  if (jet->Area()==0) return 0;
2499  //exclude jets with extremely high pt tracks which are likely misreconstructed
2500  if(jet->MaxTrackPt()>100) return 0;
2501 
2502  //passed all above cuts
2503  return 1;
2504 }
2505 
2506 //void AliAnalysisTaskEmcalJetHadEPpid::FillAnalysisSummaryHistogram() const
2508 {
2509  // fill the analysis summary histrogram, saves all relevant analysis settigns
2510  h->GetXaxis()->SetBinLabel(1, "TPC: Unidentified");
2511  h->GetXaxis()->SetBinLabel(2, "TPC: Pion");
2512  h->GetXaxis()->SetBinLabel(3, "TPC: Kaon");
2513  h->GetXaxis()->SetBinLabel(4, "TPC: Proton");
2514  h->GetXaxis()->SetBinLabel(5, "ITS: Unidentified");
2515  h->GetXaxis()->SetBinLabel(6, "ITS: Pion");
2516  h->GetXaxis()->SetBinLabel(7, "ITS: Kaon");
2517  h->GetXaxis()->SetBinLabel(8, "ITS: Proton");
2518  h->GetXaxis()->SetBinLabel(9, "TOF: Unidentified");
2519  h->GetXaxis()->SetBinLabel(10, "TOF: Pion");
2520  h->GetXaxis()->SetBinLabel(11, "TOF: Kaon");
2521  h->GetXaxis()->SetBinLabel(12, "TOF: Proton");
2522  h->GetXaxis()->SetBinLabel(14, "Unidentified tracks");
2523 
2524  // set x-axis labels vertically
2525  h->LabelsOption("v");
2526 }
2527 
2528 //void AliAnalysisTaskEmcalJetHadEPpid::FillAnalysisSummaryHistogram() const
2530 {
2531  // label bins of the analysis event summary
2532  h->GetXaxis()->SetBinLabel(1, "All events started");
2533  h->GetXaxis()->SetBinLabel(2, "object check");
2534  h->GetXaxis()->SetBinLabel(3, "aod/esd check");
2535  h->GetXaxis()->SetBinLabel(4, "centrality check");
2536  h->GetXaxis()->SetBinLabel(5, "zvertex check");
2537  h->GetXaxis()->SetBinLabel(6, "list check");
2538  h->GetXaxis()->SetBinLabel(7, "track/jet pointer check");
2539  h->GetXaxis()->SetBinLabel(8, "tracks & jets < than 1 check");
2540  h->GetXaxis()->SetBinLabel(9, "after track/jet loop to get highest pt");
2541  h->GetXaxis()->SetBinLabel(10, "accepted jets");
2542  h->GetXaxis()->SetBinLabel(11, "jets meeting pt threshold");
2543  h->GetXaxis()->SetBinLabel(12, "accepted tracks in events w/ trigger jet");
2544  h->GetXaxis()->SetBinLabel(13, "after AliVEvent & fPIDResponse");
2545  h->GetXaxis()->SetBinLabel(14, "events before event mixing");
2546  h->GetXaxis()->SetBinLabel(15, "mixed events w/ pool");
2547  h->GetXaxis()->SetBinLabel(16, "event mixing: jets");
2548  h->GetXaxis()->SetBinLabel(17, "event mixing: nMix");
2549  h->GetXaxis()->SetBinLabel(18, "event mixing: nbackground tracks");
2550  h->GetXaxis()->SetBinLabel(19, "event mixing: THE END");
2551 
2552  // set x-axis labels vertically
2553  h->LabelsOption("v");
2554 }
2555 
2557 {
2558  // label bins of the analysis trigger selection summary
2559  h->GetXaxis()->SetBinLabel(1, "no trigger");
2560  h->GetXaxis()->SetBinLabel(2, "kAny");
2561  h->GetXaxis()->SetBinLabel(3, "kAnyINT");
2562  h->GetXaxis()->SetBinLabel(4, "kMB");
2563  h->GetXaxis()->SetBinLabel(5, "kINT7");
2564  h->GetXaxis()->SetBinLabel(6, "kEMC1");
2565  h->GetXaxis()->SetBinLabel(7, "kEMC7");
2566  h->GetXaxis()->SetBinLabel(8, "kEMC8");
2567  h->GetXaxis()->SetBinLabel(9, "kEMCEJE");
2568  h->GetXaxis()->SetBinLabel(10, "kEMCEGA");
2569  h->GetXaxis()->SetBinLabel(11, "kCentral");
2570  h->GetXaxis()->SetBinLabel(12, "kSemiCentral");
2571  h->GetXaxis()->SetBinLabel(13, "kINT8");
2572  h->GetXaxis()->SetBinLabel(14, "kEMCEJE or kMB");
2573  h->GetXaxis()->SetBinLabel(15, "kEMCEGA or kMB");
2574  h->GetXaxis()->SetBinLabel(16, "kAnyINT or kMB");
2575  h->GetXaxis()->SetBinLabel(17, "kEMCEJE & (kMB or kCentral or kSemiCentral)");
2576  h->GetXaxis()->SetBinLabel(18, "kEMCEGA & (kMB or kCentral or kSemiCentral)");
2577  h->GetXaxis()->SetBinLabel(19, "kAnyINT & (kMB or kCentral or kSemiCentral)");
2578 
2579  // set x-axis labels vertically
2580  h->LabelsOption("v");
2581  //h->LabelsDeflate("X");
2582 }
2583 
2584 //______________________________________________________________________
2585 THnSparse* AliAnalysisTaskEmcalJetHadEPpid::NewTHnSparseFCorr(const char* name, UInt_t entries) {
2586  // generate new THnSparseD, axes are defined in GetDimParamsD()
2587  Int_t count = 0;
2588  UInt_t tmp = entries;
2589  while(tmp!=0){
2590  count++;
2591  tmp = tmp &~ -tmp; // clear lowest bit
2592  }
2593 
2594  TString hnTitle(name);
2595  const Int_t dim = count;
2596  Int_t nbins[dim];
2597  Double_t xmin[dim];
2598  Double_t xmax[dim];
2599 
2600  Int_t i=0;
2601  Int_t c=0;
2602  while(c<dim && i<32){
2603  if(entries&(1<<i)){
2604  TString label("");
2605  GetDimParamsCorr(i, label, nbins[c], xmin[c], xmax[c]);
2606  hnTitle += Form(";%s",label.Data());
2607  c++;
2608  }
2609 
2610  i++;
2611  }
2612  hnTitle += ";";
2613 
2614  return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
2615 } // end of NewTHnSparseF
2616 
2617 //______________________________________________________________________________________________
2618 void AliAnalysisTaskEmcalJetHadEPpid::GetDimParamsCorr(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
2619 {
2620  //stores label and binning of axis for THnSparse
2621  const Double_t pi = TMath::Pi();
2622 
2623  switch(iEntry){
2624 
2625  case 0:
2626  label = "V0 centrality (%)";
2627  nbins = 10;
2628  xmin = 0.;
2629  xmax = 100.;
2630  break;
2631 
2632  case 1:
2633  if(fcorrJetPt) { // correct jet pt
2634  label = "Jet Corrected p_{T}";
2635  nbins = 50;
2636  xmin = -50.;
2637  xmax = 200.;
2638  } else { // don't correct jet pt
2639  label = "Jet p_{T}";
2640  if(doaltPIDbinning) { // minimize bins
2641  nbins = 20;
2642  xmin = 0.;
2643  xmax = 100.;
2644  } else { // don't minimize bins
2645  nbins = 50;
2646  xmin = 0.;
2647  xmax = 250.;
2648  }
2649  }
2650  break;
2651 
2652  case 2:
2653  label = "Relative angle: Jet and Reaction Plane";
2654  nbins = 3; // (12) 48
2655  xmin = 0.;
2656  xmax = 0.5*pi;
2657  break;
2658 
2659  case 3:
2660  label = "Z-vertex";
2661  nbins = 10;
2662  xmin = -10.;
2663  xmax = 10.;
2664  break;
2665 
2666  case 4:
2667  label = "Jet p_{T} corrected with Local Rho";
2668  if(doaltPIDbinning) {
2669  nbins = 30; // 250
2670  xmin = -50.;
2671  xmax = 100.;
2672  break;
2673  } else {
2674  nbins = 50; // 250
2675  xmin = -50.;
2676  xmax = 200.;
2677  }
2678  break;
2679 
2680  case 5:
2681  label = "Jet p_{T} corrected with Global Rho";
2682  if(doaltPIDbinning) {
2683  nbins = 30; // 250
2684  xmin = -50.;
2685  xmax = 100.;
2686  break;
2687  } else {
2688  nbins = 50; // 250
2689  xmin = -50.;
2690  xmax = 200.;
2691  }
2692  break;
2693 
2694  }// end of switch
2695 } // end of Correction (ME) sparse
2696 
2697 //________________________________________________________________________
2698 //Int_t AliAnalysisTaskEmcalJetHadEPpid::AcceptFlavourJet(AliEmcalJet* fljet, Int_t NUM, Int_t NUM2, Int_t NUM3) {
2700  // Get jet if accepted for given flavour tag
2701  // If jet not accepted return 0
2702  if(!fljet) {
2703  AliError(Form("%s:Jet not found",GetName()));
2704  return 0;
2705  }
2706 
2707  Int_t flavNUM = -99;//, flavNUM2 = -99, flavNUM3 = -99; FIXME commented out to avoid compiler warning
2708  flavNUM = NUM;
2709  //flavNUM2 = NUM2;
2710  //flavNUM3 = NUM3;
2711 
2712 /*
2713  // from the AliEmcalJet class, the tagging enumerator
2714  enum EFlavourTag{
2715  kDStar = 1<<0, kD0 = 1<<1,
2716  kSig1 = 1<<2, kSig2 = 1<<3,
2717  kBckgrd1 = 1<<4, kBckgrd2 = 1<<5, kBckgrd3 = 1<<6
2718  };
2719  // bit 0 = no tag, bit 1 = Dstar, bit 2 = D0 and so forth...
2720 */
2721 
2722  // get flavour of jet (if any)
2723  Int_t flav = -999;
2724  flav = fljet->GetFlavour();
2725 
2726  // cases (for now..)
2727  // 3 = electron rich, 5 = hadron (bkgrd) rich
2728  // if flav < 1, its not tagged, so return kFALSE (0)
2729  if(flav < 1) return 0;
2730 
2731  // if flav is not equal to what we want then return kFALSE (0)
2732  //if(flav != flavNUM && flav != flavNUM2 && flav != flavNUM3) return 0;
2733  if(flav != flavNUM) return 0;
2734 
2735  // we have the flavour we want, so return kTRUE (1)
2736  //if(flav == flavNUM || flav == flavNUM2 || flav == flavNUM3) return 1;
2737  if(flav == flavNUM) return 1;
2738 
2739  // we by default have a flavour tagged jet
2740  return 1;
2741 }
2742 
2743 //________________________________________________________________________
2744 Double_t AliAnalysisTaskEmcalJetHadEPpid::EffCorrection(Double_t trackETA, Double_t trackPT, Int_t effSwitch) const {
2745  // default (current) parameters
2746 
2747 /*
2748  // min/max range of x & y (track Eta & track pt) used by Nat
2749  Double_t trETAmin = -0.9; Double_t trETAmax = 0.9;
2750  Double_t trPTmin = 0.2; Double_t trPTmax = 12.0; ??
2751 */
2752 
2753  // x-variable = track pt, y-variable = track eta
2754  Double_t x = trackPT;
2755  Double_t y = trackETA;
2756  //y = TMath::Abs(trackETA);
2757  Double_t TRefficiency = -999;
2758  Int_t runNUM = fCurrentRunNumber;
2759  Int_t runSwitchGood = -999;
2760  Int_t centbin = -99;
2761 
2762  Double_t etaaxis = 0;
2763  Double_t ptaxis = 0;
2764 
2765  if(effSwitch < 1) {
2766  // Semi-GOOD IROC C13 runlists
2767  // DO NOT USE THESE:
2768  // Runs: (169040, 169044, 169045, 169099, 169418, 169419, 169420, 169475, 169498, 169504, 169506, 169512, 169515, 169550, 169553, 169554, 169555, 169557, 169584, 169586, 169587, 169588, 169590, 169591)
2769 
2770  // Semi-GOOD OROC C08 runlists
2771  if ((runNUM == 169975 || runNUM == 169981 || runNUM == 170038 || runNUM == 170040 || runNUM == 170083 || runNUM == 170084 || runNUM == 170085 || runNUM == 170088 || runNUM == 170089 || runNUM == 170091 || runNUM == 170152 || runNUM == 170155 || runNUM == 170159 || runNUM == 170163 || runNUM == 170193 || runNUM == 170195 || runNUM == 170203 || runNUM == 170204 || runNUM == 170228 || runNUM == 170230 || runNUM == 170268 || runNUM == 170269 || runNUM == 170270 || runNUM == 170306 || runNUM == 170308 || runNUM == 170309)) runSwitchGood = 0;
2772 
2773  // GOOD runlists
2774  if ((runNUM == 167902 || runNUM == 167903 || runNUM == 167915 || runNUM == 167920 || runNUM == 167987 || runNUM == 167988 || runNUM == 168066 || runNUM == 168068 || runNUM == 168069 || runNUM == 168076 || runNUM == 168104 || runNUM == 168107 || runNUM == 168108 || runNUM == 168115 || runNUM == 168212 || runNUM == 168310 || runNUM == 168311 || runNUM == 168322 || runNUM == 168325 || runNUM == 168341 || runNUM == 168342 || runNUM == 168361 || runNUM == 168362 || runNUM == 168458 || runNUM == 168460 || runNUM == 168461 || runNUM == 168464 || runNUM == 168467 || runNUM == 168511 || runNUM == 168512 || runNUM == 168777 || runNUM == 168826 || runNUM == 168984 || runNUM == 168988 || runNUM == 168992 || runNUM == 169035 || runNUM == 169091 || runNUM == 169094 || runNUM == 169138 || runNUM == 169143 || runNUM == 169144 || runNUM == 169145 || runNUM == 169148 || runNUM == 169156 || runNUM == 169160 || runNUM == 169167 || runNUM == 169238 || runNUM == 169411 || runNUM == 169415 || runNUM == 169417 || runNUM == 169835 || runNUM == 169837 || runNUM == 169838 || runNUM == 169846 || runNUM == 169855 || runNUM == 169858 || runNUM == 169859 || runNUM == 169923 || runNUM == 169956 || runNUM == 170027 || runNUM == 170036 || runNUM == 170081)) runSwitchGood = 1;
2775 
2776  if(fCent>=0 && fCent<10) centbin = 0;
2777  else if (fCent>=10 && fCent<30) centbin = 1;
2778  else if (fCent>=30 && fCent<50) centbin = 2;
2779  else if (fCent>=50 && fCent<90) centbin = 3;
2780 
2781  if(runSwitchGood == 0 && centbin == 0) effSwitch = 2;
2782  if(runSwitchGood == 0 && centbin == 1) effSwitch = 3;
2783  if(runSwitchGood == 0 && centbin == 2) effSwitch = 4;
2784  if(runSwitchGood == 0 && centbin == 3) effSwitch = 5;
2785  if(runSwitchGood == 1 && centbin == 0) effSwitch = 6;
2786  if(runSwitchGood == 1 && centbin == 1) effSwitch = 7;
2787  if(runSwitchGood == 1 && centbin == 2) effSwitch = 8;
2788  if(runSwitchGood == 1 && centbin == 3) effSwitch = 9;
2789  }
2790 
2791  // 0-10% centrality: Semi-Good Runs
2792  Double_t p0_10SG[17] = {0.906767, 0.0754127, 1.11638, -0.0233078, 0.795454, 0.00935385, -0.000327857, 1.08903, 0.0107272, 0.443252, -0.143411, 0.965822, 0.359156, -0.581221, 1.0739, 0.00632828, 0.706356};
2793  // 10-30% centrality: Semi-Good Runs
2794  Double_t p10_30SG[17] = {0.908011, 0.0769254, 1.11912, -0.0249449, 0.741488, 0.0361252, -0.00367954, 1.10424, 0.011472, 0.452059, -0.133282, 0.980633, 0.358222, -0.620256, 1.06871, 0.00564449, 0.753168};
2795  // 30-50% centrality: Semi-Good Runs
2796  Double_t p30_50SG[17] = {0.958708, 0.0799197, 1.10817, -0.0357678, 0.75051, 0.0607808, -0.00929713, 0.998801, 0.00692244, 0.615452, -0.0480328, 0.968431, 0.321634, -0.619066, 1.03412, 0.00656201, 0.798666};
2797  // 50-90% centrality: Semi-Good Runs
2798  Double_t p50_90SG[17] = {0.944565, 0.0807258, 1.12709, -0.0324746, 0.666452, 0.0842476, -0.00963837, 1.02829, 0.00666852, 0.549625, -0.0603107, 0.981374, 0.309374, -0.619181, 1.05367, 0.005925, 0.744887};
2799 
2800  // 0-10% centrality: Good Runs
2801  Double_t p0_10G[17] = {0.971679, 0.0767571, 1.13355, -0.0274484, 0.856652, 0.00536795, 3.90795e-05, 1.06889, 0.011007, 0.447046, -0.146626, 0.919777, 0.192601, -0.268515, 1.00243, 0.00620849, 0.709477};
2802  // 10-30% centrality: Good Runs
2803  Double_t p10_30G[17] = {0.97929, 0.0776039, 1.12213, -0.0300645, 0.844722, 0.0134788, -0.0012333, 1.07955, 0.0116835, 0.456608, -0.132743, 0.930964, 0.174175, -0.267154, 0.993118, 0.00574892, 0.765256};
2804  // 30-50% centrality: Good Runs
2805  Double_t p30_50G[17] = {0.997696, 0.0816769, 1.14341, -0.0353734, 0.752151, 0.0744259, -0.0102926, 1.01561, 0.00713274, 0.57203, -0.0640248, 0.947747, 0.102007, -0.194698, 0.999164, 0.00568476, 0.7237};
2806  // 50-90% centrality: Good Runs
2807  Double_t p50_90G[17] = {0.97041, 0.0813559, 1.12151, -0.0368797, 0.709327, 0.0701501, -0.00784043, 1.06276, 0.00676173, 0.53607, -0.0703117, 0.982534, 0.0947881, -0.18073, 1.03229, 0.00580109, 0.737801};
2808 
2809  // set up a switch for different parameter values...
2810  switch(effSwitch) {
2811  case 1 :
2812  // first switch value - TRefficiency not used so = 1
2813  TRefficiency = 1.0;
2814  break;
2815 
2816  case 2 :
2817  // Parameter values for Semi-GOOD TPC (LHC11h) runs (0-10%):
2818  ptaxis = (x<2.9)*(p0_10SG[0]*exp(-pow(p0_10SG[1]/x,p0_10SG[2])) + p0_10SG[3]*x) + (x>=2.9)*(p0_10SG[4] + p0_10SG[5]*x + p0_10SG[6]*x*x);
2819  etaaxis = (y<-0.07)*(p0_10SG[7]*exp(-pow(p0_10SG[8]/TMath::Abs(y+0.91),p0_10SG[9])) + p0_10SG[10]*y) + (y>=-0.07 && y<=0.4)*(p0_10SG[11] + p0_10SG[12]*y + p0_10SG[13]*y*y) + (y>0.4)*(p0_10SG[14]*exp(-pow(p0_10SG[15]/TMath::Abs(-y+0.91),p0_10SG[16])));
2820  TRefficiency = ptaxis*etaaxis;
2821  break;
2822 
2823  case 3 :
2824  // Parameter values for Semi-GOOD TPC (LHC11h) runs (10-30%):
2825  ptaxis = (x<2.9)*(p10_30SG[0]*exp(-pow(p10_30SG[1]/x,p10_30SG[2])) + p10_30SG[3]*x) + (x>=2.9)*(p10_30SG[4] + p10_30SG[5]*x + p10_30SG[6]*x*x);
2826  etaaxis = (y<-0.07)*(p10_30SG[7]*exp(-pow(p10_30SG[8]/TMath::Abs(y+0.91),p10_30SG[9])) + p10_30SG[10]*y) + (y>=-0.07 && y<=0.4)*(p10_30SG[11] + p10_30SG[12]*y + p10_30SG[13]*y*y) + (y>0.4)*(p10_30SG[14]*exp(-pow(p10_30SG[15]/TMath::Abs(-y+0.91),p10_30SG[16])));
2827  TRefficiency = ptaxis*etaaxis;
2828  break;
2829 
2830  case 4 :
2831  // Parameter values for Semi-GOOD TPC (LHC11h) runs (30-50%):
2832  ptaxis = (x<2.9)*(p30_50SG[0]*exp(-pow(p30_50SG[1]/x,p30_50SG[2])) + p30_50SG[3]*x) + (x>=2.9)*(p30_50SG[4] + p30_50SG[5]*x + p30_50SG[6]*x*x);
2833  etaaxis = (y<-0.07)*(p30_50SG[7]*exp(-pow(p30_50SG[8]/TMath::Abs(y+0.91),p30_50SG[9])) + p30_50SG[10]*y) + (y>=-0.07 && y<=0.4)*(p30_50SG[11] + p30_50SG[12]*y + p30_50SG[13]*y*y) + (y>0.4)*(p30_50SG[14]*exp(-pow(p30_50SG[15]/TMath::Abs(-y+0.91),p30_50SG[16])));
2834  TRefficiency = ptaxis*etaaxis;
2835  break;
2836 
2837  case 5 :
2838  // Parameter values for Semi-GOOD TPC (LHC11h) runs (50-90%):
2839  ptaxis = (x<2.9)*(p50_90SG[0]*exp(-pow(p50_90SG[1]/x,p50_90SG[2])) + p50_90SG[3]*x) + (x>=2.9)*(p50_90SG[4] + p50_90SG[5]*x + p50_90SG[6]*x*x);
2840  etaaxis = (y<-0.07)*(p50_90SG[7]*exp(-pow(p50_90SG[8]/TMath::Abs(y+0.91),p50_90SG[9])) + p50_90SG[10]*y) + (y>=-0.07 && y<=0.4)*(p50_90SG[11] + p50_90SG[12]*y + p50_90SG[13]*y*y) + (y>0.4)*(p50_90SG[14]*exp(-pow(p50_90SG[15]/TMath::Abs(-y+0.91),p50_90SG[16])));
2841  TRefficiency = ptaxis*etaaxis;
2842  break;
2843 
2844  case 6 :
2845  // Parameter values for GOOD TPC (LHC11h) runs (0-10%):
2846  ptaxis = (x<2.9)*(p0_10G[0]*exp(-pow(p0_10G[1]/x,p0_10G[2])) + p0_10G[3]*x) + (x>=2.9)*(p0_10G[4] + p0_10G[5]*x + p0_10G[6]*x*x);
2847  etaaxis = (y<0.0)*(p0_10G[7]*exp(-pow(p0_10G[8]/TMath::Abs(y+0.91),p0_10G[9])) + p0_10G[10]*y) + (y>=0.0 && y<=0.4)*(p0_10G[11] + p0_10G[12]*y + p0_10G[13]*y*y) + (y>0.4)*(p0_10G[14]*exp(-pow(p0_10G[15]/TMath::Abs(-y+0.91),p0_10G[16])));
2848  TRefficiency = ptaxis*etaaxis;
2849  break;
2850 
2851  case 7 :
2852  // Parameter values for GOOD TPC (LHC11h) runs (10-30%):
2853  ptaxis = (x<2.9)*(p10_30G[0]*exp(-pow(p10_30G[1]/x,p10_30G[2])) + p10_30G[3]*x) + (x>=2.9)*(p10_30G[4] + p10_30G[5]*x + p10_30G[6]*x*x);
2854  etaaxis = (y<0.0)*(p10_30G[7]*exp(-pow(p10_30G[8]/TMath::Abs(y+0.91),p10_30G[9])) + p10_30G[10]*y) + (y>=0.0 && y<=0.4)*(p10_30G[11] + p10_30G[12]*y + p10_30G[13]*y*y) + (y>0.4)*(p10_30G[14]*exp(-pow(p10_30G[15]/TMath::Abs(-y+0.91),p10_30G[16])));
2855  TRefficiency = ptaxis*etaaxis;
2856  break;
2857 
2858  case 8 :
2859  // Parameter values for GOOD TPC (LHC11h) runs (30-50%):
2860  ptaxis = (x<2.9)*(p30_50G[0]*exp(-pow(p30_50G[1]/x,p30_50G[2])) + p30_50G[3]*x) + (x>=2.9)*(p30_50G[4] + p30_50G[5]*x + p30_50G[6]*x*x);
2861  etaaxis = (y<0.0)*(p30_50G[7]*exp(-pow(p30_50G[8]/TMath::Abs(y+0.91),p30_50G[9])) + p30_50G[10]*y) + (y>=0.0 && y<=0.4)*(p30_50G[11] + p30_50G[12]*y + p30_50G[13]*y*y) + (y>0.4)*(p30_50G[14]*exp(-pow(p30_50G[15]/TMath::Abs(-y+0.91),p30_50G[16])));
2862  TRefficiency = ptaxis*etaaxis;
2863  break;
2864 
2865  case 9 :
2866  // Parameter values for GOOD TPC (LHC11h) runs (50-90%):
2867  ptaxis = (x<2.9)*(p50_90G[0]*exp(-pow(p50_90G[1]/x,p50_90G[2])) + p50_90G[3]*x) + (x>=2.9)*(p50_90G[4] + p50_90G[5]*x + p50_90G[6]*x*x);
2868  etaaxis = (y<0.0)*(p50_90G[7]*exp(-pow(p50_90G[8]/TMath::Abs(y+0.91),p50_90G[9])) + p50_90G[10]*y) + (y>=0.0 && y<=0.4)*(p50_90G[11] + p50_90G[12]*y + p50_90G[13]*y*y) + (y>0.4)*(p50_90G[14]*exp(-pow(p50_90G[15]/TMath::Abs(-y+0.91),p50_90G[16])));
2869  TRefficiency = ptaxis*etaaxis;
2870  break;
2871 
2872  default :
2873  // no Efficiency Switch option selected.. therefore don't correct, and set eff = 1
2874  TRefficiency = 1.0;
2875  }
2876 
2877  //if(TRefficiency < 0.1) cout<<"pt: "<<x<<" eta: "<<y<<" cent: "<<fCent<<" Eff: "<<TRefficiency<<" Ptaxis: "<<ptaxis<<" etaaxis: "<<etaaxis<<endl;
2878  //if(TRefficiency > 0.90) cout<<"pt: "<<x<<" eta: "<<y<<" cent: "<<fCent<<" Eff: "<<TRefficiency<<" Ptaxis: "<<ptaxis<<" etaaxis: "<<etaaxis<<endl;
2879 
2880  return TRefficiency;
2881 }
2882 
2883 //_____________________________________________________________________________
2885 {
2886  // by default use the ep from the event header (make sure EP selection task is enabled!)
2887  Double_t a(0), b(0), c(0), d(0), e(0), f(0), g(0), h(0);
2888  vzero[0][0] = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 2, a, b);
2889  vzero[1][0] = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 2, c, d);
2890  vzero[0][1] = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 3, e, f);
2891  vzero[1][1] = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 3, g, h);
2892 }
2893 //_____________________________________________________________________________
2895 {
2896  // define some placeholders
2897  Double_t Q2[] = {-999., -999.};
2898  Double_t Q3[] = {-999., -999.};
2899 
2900  // for all other types use calibrated event plane from the event header
2901  //
2902  // note that the code is a bit messy here, for 11h done all in this function
2903  // reason is that the procedure is much shorter as the calibration is done in another task
2904  //
2905  // define some pleaceholders to the values by reference
2906  Double_t qx2a(0), qy2a(0), qx2c(0), qy2c(0), qx3a(0), qy3a(0), qx3c(0), qy3c(0);
2907  // get the q-vectors by reference
2908  InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 2, qx2a, qy2a); // 2nd order VZERO Aside
2909  InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 2, qx2c, qy2c); // 2nd order VZERO Cside
2910  InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 3, qx3a, qy3a); // 3rd order VZERO Aside
2911  InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 3, qx3c, qy3c); // 3rd order VZERO Cside
2912 
2913  // get cache index and retrieve the chi weights for this centrality
2914 // Int_t VZEROcentralityBin(GetCentBin(fCent)); // won't work because I use different bins then calibration value bins
2915  Int_t VZEROcentralityBin(GetVZEROCentralityBin());
2916  Double_t chi2A(fChi2A->At(VZEROcentralityBin));
2917  Double_t chi2C(fChi2C->At(VZEROcentralityBin));
2918  Double_t chi3A(fChi3A->At(VZEROcentralityBin));
2919  Double_t chi3C(fChi3C->At(VZEROcentralityBin));
2920 
2921  // combine the vzeroa and vzeroc signal
2922  Q2[0] = chi2A*chi2A*qx2a+chi2C*chi2C*qx2c;
2923  Q2[1] = chi2A*chi2A*qy2a+chi2C*chi2C*qy2c;
2924  Q3[0] = chi3A*chi3A*qx3a+chi3C*chi3C*qx3c;
2925  Q3[1] = chi3A*chi3A*qy3a+chi3C*chi3C*qy3c;
2926 
2927  comb[0] = .5*TMath::ATan2(Q2[1], Q2[0]);
2928  comb[1] = (1./3.)*TMath::ATan2(Q3[1], Q3[0]);
2929 }
2930 
2931 //_____________________________________________________________________________
2933 {
2934  // grab the TPC event plane
2935  fNAcceptedTracks = 0; // reset the track counter
2936  Double_t qx2(0), qy2(0); // for psi2
2937  Double_t qx3(0), qy3(0); // for psi3
2938 
2939  // first check for track object
2940  if(fTracks) {
2941  Float_t excludeInEta = -999;
2942  if(fExcludeLeadingJetsFromFit > 0 ) { // remove the leading jet from ep estimate
2943  if(fLeadingJet) excludeInEta = fLeadingJet->Eta();
2944  }
2945 
2946  // loop over tracks in TPC
2947  for(Int_t iTPC(0); iTPC < fTracks->GetEntries(); iTPC++) {
2948  AliVParticle *track = dynamic_cast<AliVParticle*>(fTracks->At(iTPC));
2949 
2950  // apply general track cuts
2951  if(TMath::Abs(track->Eta())>0.9) continue;
2952  if(track->Pt()<0.15) continue;
2953  if(track->Pt() < fSoftTrackMinPt_ep || track->Pt() > fSoftTrackMaxPt_ep) continue; // APPLY soft cut
2954  if(fExcludeLeadingJetsFromFit > 0 &&( (TMath::Abs(track->Eta() - excludeInEta) < fJetRad*fExcludeLeadingJetsFromFit ) || (TMath::Abs(track->Eta()) - fJetRad - fEtamax ) > 0 )) continue;
2955  fNAcceptedTracks++;
2956 
2957  // sum up q-vectors
2958  qx2+= TMath::Cos(2.*track->Phi());
2959  qy2+= TMath::Sin(2.*track->Phi());
2960  qx3+= TMath::Cos(3.*track->Phi());
2961  qy3+= TMath::Sin(3.*track->Phi());
2962  }
2963  }
2964  tpc[0] = .5*TMath::ATan2(qy2, qx2);
2965  tpc[1] = (1./3.)*TMath::ATan2(qy3, qx3);
2966 }
2967 
2968 //_____________________________________________________________________________
2969 void AliAnalysisTaskEmcalJetHadEPpid::CalculateEventPlaneResolution(Double_t vzero[2][2], Double_t* vzeroComb, Double_t* tpc)
2970 {
2971  // fill the profiles for the resolution parameters
2972  Int_t centbin = GetCentBin(fCent);
2973  fInCentralitySelection = centbin;
2974 
2975  // R2 resolution for 2nd order event plane
2976  fProfV2Resolution[fInCentralitySelection]->Fill(2., TMath::Cos(2.*(vzero[0][0] - vzero[1][0])));
2977  fProfV2Resolution[fInCentralitySelection]->Fill(3., TMath::Cos(2.*(vzero[1][0] - vzero[0][0])));
2978  fProfV2Resolution[fInCentralitySelection]->Fill(4., TMath::Cos(2.*(vzero[0][0] - tpc[0])));
2979  fProfV2Resolution[fInCentralitySelection]->Fill(5., TMath::Cos(2.*(tpc[0] - vzero[0][0])));
2980  fProfV2Resolution[fInCentralitySelection]->Fill(6., TMath::Cos(2.*(vzero[1][0] - tpc[0])));
2981  fProfV2Resolution[fInCentralitySelection]->Fill(7., TMath::Cos(2.*(tpc[0] - vzero[1][0])));
2982  // R3 resolution for 2nd order event plane
2983  fProfV3Resolution[fInCentralitySelection]->Fill(2., TMath::Cos(3.*(vzero[0][0] - vzero[1][0])));
2984  fProfV3Resolution[fInCentralitySelection]->Fill(3., TMath::Cos(3.*(vzero[1][0] - vzero[0][0])));
2985  fProfV3Resolution[fInCentralitySelection]->Fill(4., TMath::Cos(3.*(vzero[0][0] - tpc[0])));
2986  fProfV3Resolution[fInCentralitySelection]->Fill(5., TMath::Cos(3.*(tpc[0] - vzero[0][0])));
2987  fProfV3Resolution[fInCentralitySelection]->Fill(6., TMath::Cos(3.*(vzero[1][0] - tpc[0])));
2988  fProfV3Resolution[fInCentralitySelection]->Fill(7., TMath::Cos(3.*(tpc[0] - vzero[1][0])));
2989  // R4 resolution for 2nd order event plane
2990  fProfV4Resolution[fInCentralitySelection]->Fill(2., TMath::Cos(4.*(vzero[0][0] - vzero[1][0])));
2991  fProfV4Resolution[fInCentralitySelection]->Fill(3., TMath::Cos(4.*(vzero[1][0] - vzero[0][0])));
2992  fProfV4Resolution[fInCentralitySelection]->Fill(4., TMath::Cos(4.*(vzero[0][0] - tpc[0])));
2993  fProfV4Resolution[fInCentralitySelection]->Fill(5., TMath::Cos(4.*(tpc[0] - vzero[0][0])));
2994  fProfV4Resolution[fInCentralitySelection]->Fill(6., TMath::Cos(4.*(vzero[1][0] - tpc[0])));
2995  fProfV4Resolution[fInCentralitySelection]->Fill(7., TMath::Cos(4.*(tpc[0] - vzero[1][0])));
2996  // R5 resolution for 2nd order event plane
2997  fProfV5Resolution[fInCentralitySelection]->Fill(2., TMath::Cos(5.*(vzero[0][0] - vzero[1][0])));
2998  fProfV5Resolution[fInCentralitySelection]->Fill(3., TMath::Cos(5.*(vzero[1][0] - vzero[0][0])));
2999  fProfV5Resolution[fInCentralitySelection]->Fill(4., TMath::Cos(5.*(vzero[0][0] - tpc[0])));
3000  fProfV5Resolution[fInCentralitySelection]->Fill(5., TMath::Cos(5.*(tpc[0] - vzero[0][0])));
3001  fProfV5Resolution[fInCentralitySelection]->Fill(6., TMath::Cos(5.*(vzero[1][0] - tpc[0])));
3002  fProfV5Resolution[fInCentralitySelection]->Fill(7., TMath::Cos(5.*(tpc[0] - vzero[1][0])));
3003 
3004  // for the resolution of the combined vzero event plane, use two tpc halves as uncorrelated subdetectors
3005  Double_t qx2a(0), qy2a(0); // for psi2a, negative eta
3006  Double_t qx3a(0), qy3a(0); // for psi3a, negative eta
3007  Double_t qx2b(0), qy2b(0); // for psi2c, positive eta
3008  Double_t qx3b(0), qy3b(0); // for psi3c, positive eta
3009  // not needed
3010  Double_t qx4a(0), qy4a(0); // for psi4a, negative eta
3011  Double_t qx5a(0), qy5a(0); // for psi4a, negative eta
3012  Double_t qx4b(0), qy4b(0); // for psi5c, positive eta
3013  Double_t qx5b(0), qy5b(0); // for psi5c, positive eta
3014 
3015  // check for track object and loop over tpc tracks
3016  if(fTracks) {
3017  Int_t iTracks(fTracks->GetEntriesFast());
3018  for(Int_t iTPC(0); iTPC < iTracks; iTPC++) {
3019  AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(iTPC));
3020 
3021  // apply track cuts
3022  if(TMath::Abs(track->Eta())>0.9) continue;
3023  if(track->Pt()<0.15) continue;
3024  if(track->Pt() < fSoftTrackMinPt_ep || track->Pt() > fSoftTrackMaxPt_ep) continue;
3025  if(track->Eta() < 0 ) {
3026  // negative Eta q-vect calculation
3027  qx2a+= TMath::Cos(2.*track->Phi());
3028  qy2a+= TMath::Sin(2.*track->Phi());
3029  qx3a+= TMath::Cos(3.*track->Phi());
3030  qy3a+= TMath::Sin(3.*track->Phi());
3031  qx4a+= TMath::Cos(4.*track->Phi());
3032  qy4a+= TMath::Sin(4.*track->Phi());
3033  qx5a+= TMath::Cos(5.*track->Phi());
3034  qy5a+= TMath::Sin(5.*track->Phi());
3035  } else if (track->Eta() > 0) {
3036  // positive Eta q-vect calculation
3037  qx2b+= TMath::Cos(2.*track->Phi());
3038  qy2b+= TMath::Sin(2.*track->Phi());
3039  qx3b+= TMath::Cos(3.*track->Phi());
3040  qy3b+= TMath::Sin(3.*track->Phi());
3041  qx4b+= TMath::Cos(4.*track->Phi());
3042  qy4b+= TMath::Sin(4.*track->Phi());
3043  qx5b+= TMath::Cos(5.*track->Phi());
3044  qy5b+= TMath::Sin(5.*track->Phi());
3045  }
3046  }
3047  }
3048 
3049  // combine x-y vectors for neg and pos Eta ranges
3050  Double_t tpca2(.5*TMath::ATan2(qy2a, qx2a));
3051  Double_t tpcb2(.5*TMath::ATan2(qy2b, qx2b));
3052 
3053  // not used, but can be for 3rd, 4th, and 5th order event planes
3054 /*
3055  Double_t tpca3((1./3.)*TMath::ATan2(qy3a, qx3a));
3056  Double_t tpcb3((1./3.)*TMath::ATan2(qy3b, qx3b));
3057  Double_t tpca4(.25*TMath::ATan2(qy4a, qx4a));
3058  Double_t tpcb4(.25*TMath::ATan2(qy4b, qx4b));
3059  Double_t tpca5((.2)*TMath::ATan2(qy5a, qx5a));
3060  Double_t tpcb5((.2)*TMath::ATan2(qy5b, qx5b));
3061 */
3062 
3063  // R2 resolution for 2nd order event plane
3064  fProfV2Resolution[fInCentralitySelection]->Fill(8., TMath::Cos(2.*(vzeroComb[0] - tpca2)));
3065  fProfV2Resolution[fInCentralitySelection]->Fill(9., TMath::Cos(2.*(vzeroComb[0] - tpcb2)));
3066  fProfV2Resolution[fInCentralitySelection]->Fill(10., TMath::Cos(2.*(tpca2 - tpcb2)));
3067  // R3 resolution for 2nd order event plane
3068  fProfV3Resolution[fInCentralitySelection]->Fill(8., TMath::Cos(3.*(vzeroComb[1] - tpca2)));
3069  fProfV3Resolution[fInCentralitySelection]->Fill(9., TMath::Cos(3.*(vzeroComb[1] - tpcb2)));
3070  fProfV3Resolution[fInCentralitySelection]->Fill(10., TMath::Cos(3.*(tpca2 - tpcb2)));
3071  // R4 resolution for 2nd order event plane
3072  fProfV4Resolution[fInCentralitySelection]->Fill(8., TMath::Cos(4.*(vzeroComb[0] - tpca2)));
3073  fProfV4Resolution[fInCentralitySelection]->Fill(9., TMath::Cos(4.*(vzeroComb[0] - tpcb2)));
3074  fProfV4Resolution[fInCentralitySelection]->Fill(10., TMath::Cos(4.*(tpca2 - tpcb2)));
3075  // R5 resolution for 2nd order event plane
3076  fProfV5Resolution[fInCentralitySelection]->Fill(8., TMath::Cos(5.*(vzeroComb[1] - tpca2)));
3077  fProfV5Resolution[fInCentralitySelection]->Fill(9., TMath::Cos(5.*(vzeroComb[1] - tpcb2)));
3078  fProfV5Resolution[fInCentralitySelection]->Fill(10., TMath::Cos(5.*(tpca2 - tpcb2)));
3079 }
3080 
3081 //_____________________________________________________________________________i
3083 {
3084  // make sure calibration parameters are present for vzero ep 11h data
3085  //
3086  // chi values can be calculated using the static helper function
3087  // AliAnalysisTaskJetV2::CalculateEventPlaneChi(Double_t res) where res is the event plane
3088  // resolution in a given centrality bin
3089  //
3090  // the resolutions that were used for these defaults are
3091  // this might need a bit of updating as they were read 'by-eye' from a performance plot ..
3092  // Double_t R2VZEROA[] = {.35, .40, .48, .50, .48, .45, .38, .26, .16};
3093  // Double_t R2VZEROC[] = {.45, .60, .70, .73, .68, .60, .40, .36, .17};
3094  // Double_t R3VZEROA[] = {.22, .23, .22, .19, .15, .12, .08, .00, .00};
3095  // Double_t R3VZEROC[] = {.30, .30, .28, .25, .22, .17, .11, .00, .00};
3096 
3097  Double_t chiA2[] = {0.582214, 0.674622, 0.832214, 0.873962, 0.832214, 0.771423, 0.637146, 0.424255, 0.257385};
3098  Double_t chiC2[] = {0.771423, 1.10236, 1.38116, 1.48077, 1.31964, 1.10236, 0.674622, 0.600403, 0.273865};
3099  Double_t chiA3[] = {0.356628, 0.373474, 0.356628, 0.306702, 0.24115, 0.192322, 0.127869, 6.10352e-05, 6.10352e-05};
3100  Double_t chiC3[] = {0.493347, 0.493347, 0.458557, 0.407166, 0.356628, 0.273865, 0.176208, 6.10352e-05, 6.10352e-05};
3101 
3102  // when the user wants to, set the weights to 1 (effectively disabling them)
3103  if(!fUseChiWeightForVZERO) {
3104  for(Int_t i(0); i < 9; i++) {
3105  chiA2[i] = 1.;
3106  chiA3[i] = 1.;
3107  chiC2[i] = 1.;
3108  chiC3[i] = 1.;
3109  }
3110  }
3111 
3112  if(!fChi2A) fChi2A = new TArrayD(9, chiA2);
3113  if(!fChi2C) fChi2C = new TArrayD(9, chiC2);
3114  if(!fChi3A) fChi3A = new TArrayD(9, chiA3);
3115  if(!fChi3C) fChi3C = new TArrayD(9, chiC3);
3116 }
3117 
3118 //_____________________________________________________________________________
3120 {
3121  // return cache index number corresponding to the event centrality
3122  Float_t v0Centr(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"));
3123  if(v0Centr < 5) return 0;
3124  else if(v0Centr < 10) return 1;
3125  else if(v0Centr < 20) return 2;
3126  else if(v0Centr < 30) return 3;
3127  else if(v0Centr < 40) return 4;
3128  else if(v0Centr < 50) return 5;
3129  else if(v0Centr < 60) return 6;
3130  else if(v0Centr < 70) return 7;
3131  else return 8;
3132 }
3133 
3134 //_____________________________________________________________________________
3136  // return pointer to the highest pt jet (before background subtraction) within acceptance
3137  // only rudimentary cuts are applied on this level, hence the implementation outside of
3138  // the framework
3139  if(fJets) {
3140  Int_t iJets(fJets->GetEntriesFast());
3141  Double_t pt(0);
3142  AliEmcalJet* leadingJet(0x0);
3143  if(!localRho) {
3144  for(Int_t i(0); i < iJets; i++) {
3145  AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(i));
3146  if(!AcceptMyJet(jet)) continue;
3147  if(jet->Pt() > pt) {
3148  leadingJet = jet;
3149  pt = leadingJet->Pt();
3150  }
3151  }
3152  return leadingJet;
3153  } else {
3154  // return leading jet after background subtraction
3155  Double_t rho(0);
3156  for(Int_t i(0); i < iJets; i++) {
3157  AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(i));
3158  if(!AcceptMyJet(jet)) continue;
3159  //rho = localRho->GetLocalVal(jet->Phi(), GetJetContainer()->GetJetRadius(), localRho->GetVal());
3160  rho = localRho->GetLocalVal(jet->Phi(), fJetRad);
3161  if((jet->Pt()-jet->Area()*rho) > pt) {
3162  leadingJet = jet;
3163  pt = (leadingJet->Pt()-jet->Area()*rho);
3164  }
3165  }
3166  return leadingJet;
3167 
3168  }
3169  }
3170 
3171  return 0x0;
3172 }
3173 
3174 //_____________________________________________________________________________
3176 {
3177  // return chi for given resolution to combine event plane estimates from two subevents
3178  // see Phys. Rev. C no. CS6346 (http://arxiv.org/abs/nucl-ex/9805001)
3179  Double_t chi(2.), delta(1.), con((TMath::Sqrt(TMath::Pi()))/(2.*TMath::Sqrt(2)));
3180  for (Int_t i(0); i < 15; i++) {
3181  chi = ((con*chi*TMath::Exp(-chi*chi/4.)*(TMath::BesselI0(chi*chi/4.)+TMath::BesselI1(chi*chi/4.))) < res) ? chi + delta : chi - delta;
3182  delta = delta / 2.;
3183  }
3184  return chi;
3185 }
3186 
3188  // check and fill a Event Selection QA histogram for different trigger selections after cuts
3189  if(trig == 0) h->Fill(1);
3190  if(trig & AliVEvent::kAny) h->Fill(2);
3191  if(trig & AliVEvent::kAnyINT) h->Fill(3);
3192  if(trig & AliVEvent::kMB) h->Fill(4);
3193  if(trig & AliVEvent::kINT7) h->Fill(5);
3194  if(trig & AliVEvent::kEMC1) h->Fill(6);
3195  if(trig & AliVEvent::kEMC7) h->Fill(7);
3196  if(trig & AliVEvent::kEMC8) h->Fill(8);
3197  if(trig & AliVEvent::kEMCEJE) h->Fill(9);
3198  if(trig & AliVEvent::kEMCEGA) h->Fill(10);
3199  if(trig & AliVEvent::kCentral) h->Fill(11);
3200  if(trig & AliVEvent::kSemiCentral) h->Fill(12);
3201  if(trig & AliVEvent::kINT8) h->Fill(13);
3202 
3203  if(trig & (AliVEvent::kEMCEJE | AliVEvent::kMB)) h->Fill(14);
3204  if(trig & (AliVEvent::kEMCEGA | AliVEvent::kMB)) h->Fill(15);
3205  if(trig & (AliVEvent::kAnyINT | AliVEvent::kMB)) h->Fill(16);
3206 
3207  if(trig & (AliVEvent::kEMCEJE & (AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral))) h->Fill(17);
3208  if(trig & (AliVEvent::kEMCEGA & (AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral))) h->Fill(18);
3209  if(trig & (AliVEvent::kAnyINT & (AliVEvent::kMB | AliVEvent::kCentral | AliVEvent::kSemiCentral))) h->Fill(19);
3210 
3211  // label bins of the analysis trigger selection summary
3212  h->GetXaxis()->SetBinLabel(1, "no trigger");
3213  h->GetXaxis()->SetBinLabel(2, "kAny");
3214  h->GetXaxis()->SetBinLabel(3, "kAnyINT");
3215  h->GetXaxis()->SetBinLabel(4, "kMB");
3216  h->GetXaxis()->SetBinLabel(5, "kINT7");
3217  h->GetXaxis()->SetBinLabel(6, "kEMC1");
3218  h->GetXaxis()->SetBinLabel(7, "kEMC7");
3219  h->GetXaxis()->SetBinLabel(8, "kEMC8");
3220  h->GetXaxis()->SetBinLabel(9, "kEMCEJE");
3221  h->GetXaxis()->SetBinLabel(10, "kEMCEGA");
3222  h->GetXaxis()->SetBinLabel(11, "kCentral");
3223  h->GetXaxis()->SetBinLabel(12, "kSemiCentral");
3224  h->GetXaxis()->SetBinLabel(13, "kINT8");
3225  h->GetXaxis()->SetBinLabel(14, "kEMCEJE or kMB");
3226  h->GetXaxis()->SetBinLabel(15, "kEMCEGA or kMB");
3227  h->GetXaxis()->SetBinLabel(16, "kAnyINT or kMB");
3228  h->GetXaxis()->SetBinLabel(17, "kEMCEJE & (kMB or kCentral or kSemiCentral)");
3229  h->GetXaxis()->SetBinLabel(18, "kEMCEGA & (kMB or kCentral or kSemiCentral)");
3230  h->GetXaxis()->SetBinLabel(19, "kAnyINT & (kMB or kCentral or kSemiCentral)");
3231 
3232  // set x-axis labels vertically
3233  h->LabelsOption("v");
3234  //h->LabelsDeflate("X");
3235 
3236  return h;
3237 }
TObjArray * CloneAndReduceTrackList(TObjArray *tracks)
Int_t charge
Double_t Area() const
Definition: AliEmcalJet.h:69
virtual AliVParticle * GetNextAcceptParticle()
virtual Int_t GetpTtrackBin(Double_t pt) const
Double_t RelativeEPJET(Double_t jetAng, Double_t EPAng) const
TProfile * fProfV5Resolution[10]
resolution parameters for v4
AliEmcalJet * GetLeadingJet(AliLocalRhoParameter *localRho=0x0)
THnSparse * fhnJH
// mixed events matrix
AliEmcalJet * fLeadingJet
number of accepted tracks
ClassImp(AliAnalysisTaskEmcalJetHadEPpid) AliAnalysisTaskEmcalJetHadEPpid
Double_t Eta() const
Definition: AliEmcalJet.h:60
Int_t GetFlavour() const
Definition: AliEmcalJet.h:259
static Double_t CalculateEventPlaneChi(Double_t res)
TArrayD * fChi2A
resolution parameters for v5
Double_t Phi() const
Definition: AliEmcalJet.h:55
Double_t GetLocalVal(Double_t phi, Double_t r, Double_t n) const
Double_t fEPV0
!event plane V0
TList * list
TH1 * fHistTrackPtallcent
// phi distrubtion of mixed events
Double_t GetLeadingHadronPt(AliEmcalJet *jet, Int_t c=0)
virtual Int_t GetEtaBin(Double_t eta) const
TProfile * fProfV4Resolution[10]
resolution parameters for v3
AliLocalRhoParameter * GetLocalRhoFromEvent(const char *name)
void CalculateEventPlaneResolution(Double_t vzero[2][2], Double_t *vzeroComb, Double_t *tpc)
virtual void GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
TList * fOutput
!output list
virtual Int_t GetpTjetBin(Double_t pt) const
Double_t fEPV0C
!event plane V0C
UShort_t GetNumberOfConstituents() const
Definition: AliEmcalJet.h:84
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:83
AliParticleContainer * GetParticleContainer(Int_t i=0) const
virtual void GetDimParamsPID(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
Double_t Eta() const
Definition: AliPicoTrack.h:37
AliRhoParameter * GetRhoFromEvent(const char *name)
virtual THnSparse * NewTHnSparseFCorr(const char *name, UInt_t entries)
UShort_t GetNumberOfClusters() const
Definition: AliEmcalJet.h:82
TH1F * fHistJetHaddPHI
// mixed events phi-eta distributions
TProfile * fProfV3Resolution[10]
resolution parameters for v2
virtual void GetDimParamsCorr(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
TH2F * fHistJetPtvsTrackPt[6]
number of jets versus Centrality
const Double_t ptmax
AliRhoParameter * fRho
!event rho
virtual THnSparse * NewTHnSparseF(const char *name, UInt_t entries)
Double_t MaxTrackPt() const
Definition: AliEmcalJet.h:99
Double_t Phi() const
Definition: AliPicoTrack.h:33
AliPIDResponse * fPIDResponse
// event pool Manager object
Double_t fCent
!event centrality
TClonesArray * GetArray() const
virtual THnSparse * NewTHnSparseFPID(const char *name, UInt_t entries)
AliLocalRhoParameter * fLocalRho
!local event rho
Double_t Pt() const
Definition: AliPicoTrack.h:23
AliClusterContainer * fCaloClustersCont
Tracks.
Double_t RelativePhi(Double_t mphi, Double_t vphi) const
virtual Int_t GetzVertexBin(Double_t zVtx) const
void SetClassName(const char *clname)
TClonesArray * fJets
jet collection array
Double_t Pt() const
Definition: AliEmcalJet.h:47
Double_t EffCorrection(Double_t trkETA, Double_t trkPT, Int_t effswitch) const
const char * GetName() const
virtual Int_t AcceptFlavourJet(AliEmcalJet *jet, Int_t NUM)
TClonesArray * fTracks
!tracks
TProfile * fProfV2Resolution[10]
// sparse to get # jet triggers
Bool_t fCreateHisto
whether or not create histograms
void CalculateEventPlaneCombinedVZERO(Double_t *comb) const
void SetMakeGeneralHistograms(Bool_t g)
Double_t fEPV0A
!event plane V0A
Double_t MaxClusterPt() const
Definition: AliEmcalJet.h:98
const Double_t pi
const Int_t nbins
THnSparse * fhnCorr
// jet hadron events matrix
TH2F * fHistMEphieta
// single events phi-eta distributions
Double_t fRhoVal
!event rho value, same for local rho
void ResetCurrentID(Int_t i=-1)
virtual Int_t GetCentBin(Double_t cent) const
void CalculateEventPlaneVZERO(Double_t vzero[2][2]) const
Short_t Charge() const
Definition: AliPicoTrack.h:39