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