AliPhysics  fe039ad (fe039ad)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskChargedJetsHadronCF.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2016, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: R. Haake. *
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 #include <algorithm>
17 #include <vector>
18 #include <TClonesArray.h>
19 #include <TF1.h>
20 #include <TH1F.h>
21 #include <TH2F.h>
22 #include <TH3F.h>
23 #include <THn.h>
24 #include <TTree.h>
25 #include <TList.h>
26 #include <TLorentzVector.h>
27 
28 #include "AliEmcalPythiaInfo.h"
29 #include "AliMCEvent.h"
30 #include "AliPythia.h"
31 #include "AliStack.h"
32 
33 #include "AliVTrack.h"
34 #include "AliVHeader.h"
35 #include "AliEmcalJet.h"
36 #include "AliRhoParameter.h"
37 #include "AliLog.h"
38 #include "AliJetContainer.h"
39 #include "AliTrackContainer.h"
40 #include "AliAODTrack.h"
41 #include "AliAODMCParticle.h"
42 #include "AliAODPid.h"
43 #include "AliPicoTrack.h"
44 #include "AliVParticle.h"
45 #include "TRandom3.h"
47 #include "AliBasicParticle.h"
48 
50 
51 //________________________________________________________________________
53 {
54 // dummy destructor
55 }
56 
60 //________________________________________________________________________
63  fJetsCont(0),
64  fTracksCont(0),
65  fEventExtractionPercentage(0),
66  fEventExtractionMinJetPt(0),
67  fEventExtractionMaxJetPt(0),
68  fConstPtFilterBit(1024),
69  fNumberOfCentralityBins(10),
70  fJetsOutput(),
71  fTracksOutput(),
72  fJetParticleArrayName("JetsDPhiBasicParticles"),
73  fTrackParticleArrayName(""),
74  fJetEmbeddingArray(),
75  fJetEmbeddingArrayName(""),
76  fJetEmbeddingTrackArrayName(""),
77  fJetEmbeddingMaxDistance(0.3),
78  fJetEmbeddingNumMatchedJets(2),
79  fJetEmbeddingUsePerTrackMCPercentage(kTRUE),
80  fJetEmbeddingUseBgrdForMCPercentage(kFALSE),
81  fJetEmbeddingCreatePtPlotPerCut(kFALSE),
82  fJetEmbeddingCuts(),
83  fJetVetoArray(),
84  fJetVetoArrayName(""),
85  fJetVetoJetByJet(1),
86  fMatchedJets(),
87  fRandom(0),
88  fTracksTree(0),
89  fTreeBufferTrack(0),
90  fTreeBufferPID(0),
91  fTreeBufferPDG(0),
92  fTrackExtractionPercentagePower(0),
93  fNumRandomConesPerEvent(10),
94  fUseDataConstituents(kTRUE),
95  fUseMCConstituents(kTRUE),
96  fJetOutputMode(0),
97  fLeadingJet(0),
98  fSubleadingJet(0),
99  fInitialPartonMatchedJet1(0),
100  fInitialPartonMatchedJet2(0),
101  fAcceptedJets(0),
102  fAcceptedTracks(0)
103 {
104  // Default constructor.
105  SetMakeGeneralHistograms(kTRUE);
106  fRandom = new TRandom3(0);
107 }
108 
109 
110 //________________________________________________________________________
112  AliAnalysisTaskEmcalJet(name, kTRUE),
113  fJetsCont(0),
114  fTracksCont(0),
115  fEventExtractionPercentage(0),
116  fEventExtractionMinJetPt(0),
117  fEventExtractionMaxJetPt(0),
118  fConstPtFilterBit(1024),
119  fNumberOfCentralityBins(10),
120  fJetsOutput(),
121  fTracksOutput(),
122  fJetParticleArrayName("JetsDPhiBasicParticles"),
123  fTrackParticleArrayName(""),
124  fJetEmbeddingArray(),
125  fJetEmbeddingArrayName(""),
126  fJetEmbeddingTrackArrayName(""),
127  fJetEmbeddingMaxDistance(0.3),
128  fJetEmbeddingNumMatchedJets(2),
129  fJetEmbeddingUsePerTrackMCPercentage(kTRUE),
130  fJetEmbeddingUseBgrdForMCPercentage(kFALSE),
131  fJetEmbeddingCreatePtPlotPerCut(kFALSE),
132  fJetEmbeddingCuts(),
133  fJetVetoArray(),
134  fJetVetoArrayName(""),
135  fJetVetoJetByJet(1),
136  fMatchedJets(),
137  fRandom(0),
138  fTracksTree(0),
139  fTreeBufferTrack(0),
140  fTreeBufferPID(0),
141  fTreeBufferPDG(0),
142  fTrackExtractionPercentagePower(0),
143  fNumRandomConesPerEvent(10),
144  fUseDataConstituents(kTRUE),
145  fUseMCConstituents(kTRUE),
146  fJetOutputMode(0),
147  fLeadingJet(0),
148  fSubleadingJet(0),
149  fInitialPartonMatchedJet1(0),
150  fInitialPartonMatchedJet2(0),
151  fAcceptedJets(0),
152  fAcceptedTracks(0)
153 {
154  // Constructor
156  fRandom = new TRandom3(0);
157 }
158 
159 //________________________________________________________________________
161 {
162  // Destructor.
163 }
164 
165 //________________________________________________________________________
167 {
169 
170  // ### Basic container settings
172  if(fJetsCont) { //get particles connected to jets
173  fJetsCont->PrintCuts();
175  } else { //no jets, just analysis tracks
177  }
178  if(fTracksCont) fTracksCont->SetClassName("AliVTrack");
179 
180  // ### Create all histograms
181  fHistEventRejection->GetXaxis()->SetBinLabel(15,"JetVeto");
182 
183  // Track QA plots
184  AddHistogram2D<TH2D>("hTrackPt", "Tracks p_{T} distribution", "", 300, 0., 300., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)", "Centrality", "dN^{Tracks}/dp_{T}");
185  AddHistogram2D<TH2D>("hTrackPhi", "Track angular distribution in #phi", "LEGO2", 180, 0., 2*TMath::Pi(), fNumberOfCentralityBins, 0, 100, "#phi", "Centrality", "dN^{Tracks}/(d#phi)");
186  AddHistogram2D<TH2D>("hTrackEta", "Track angular distribution in #eta", "LEGO2", 100, -2.5, 2.5, fNumberOfCentralityBins, 0, 100, "#eta", "Centrality", "dN^{Tracks}/(d#eta)");
187  AddHistogram2D<TH2D>("hTrackPhiEta", "Track angular distribution #phi/#eta", "COLZ", 180, 0., 2*TMath::Pi(), 100, -2.5, 2.5, "#phi", "#eta", "dN^{Tracks}/d#phi d#eta");
188 
189  AddHistogram2D<TH2D>("hLeadingTrackPt", "Leading tracks p_{T} distribution", "", 300, 0., 300., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)", "Centrality", "dN^{Tracks}/dp_{T}");
190  AddHistogram2D<TH2D>("hLeadingTrackPhi", "Leading tracks angular distribution in #phi", "LEGO2", 180, 0., 2*TMath::Pi(), fNumberOfCentralityBins, 0, 100, "#phi", "Centrality", "dN^{Tracks}/(d#phi)");
191  AddHistogram2D<TH2D>("hLeadingTrackEta", "Leading tracks angular distribution in #eta", "LEGO2", 100, -2.5, 2.5, fNumberOfCentralityBins, 0, 100, "#eta", "Centrality", "dN^{Tracks}/(d#eta)");
192  AddHistogram2D<TH2D>("hLeadingTrackPhiEta", "Track angular distribution #phi/#eta", "COLZ", 180, 0., 2*TMath::Pi(), 100, -2.5, 2.5, "#phi", "#eta", "dN^{Tracks}/d#phi d#eta");
193 
194  AddHistogram2D<TH2D>("hTrackEtaPt", "Track angular distribution in #eta vs. p_{T}", "LEGO2", 100, -2.5, 2.5, 300, 0., 300., "#eta", "p_{T} (GeV/c)", "dN^{Tracks}/(d#eta dp_{T})");
195  AddHistogram2D<TH2D>("hTrackPhiPt", "Track angular distribution in #phi vs. p_{T}", "LEGO2", 180, 0, 2*TMath::Pi(), 300, 0., 300., "#phi", "p_{T} (GeV/c)", "dN^{Tracks}/(d#phi dp_{T})");
196 
197  // Create plots for each embedding cut
198  for(Int_t i = -1; i<static_cast<Int_t>(fJetEmbeddingCuts.size()); i++)
199  {
200  const char* appendix = "";
201  if(i>-1)
202  {
204  appendix = Form("_%s", currentCut.fCutName.Data());
205 
206  // Don't double-add cuts
207  if( static_cast<TH1*>(fOutput->FindObject(Form("hJetPtRaw%s", appendix))) )
208  continue;
209  }
210  // Jet QA plots
211  AddHistogram2D<TH2D>(Form("hJetPtRaw%s", appendix), "Jets p_{T} distribution (no bgrd. corr.)", "", 300, 0., 300., fNumberOfCentralityBins, 0, 100, "p_{T, jet} (GeV/c)", "Centrality", "dN^{Jets}/dp_{T}");
212  AddHistogram2D<TH2D>(Form("hJetPt%s", appendix), "Jets p_{T} distribution (background subtracted)", "", 400, -100., 300., fNumberOfCentralityBins, 0, 100, "p_{T, jet} (GeV/c)", "Centrality", "dN^{Jets}/dp_{T}");
213  AddHistogram2D<TH2D>(Form("hJetPhi%s", appendix), "Jet angular distribution #phi", "LEGO2", 180, 0., 2*TMath::Pi(), fNumberOfCentralityBins, 0, 100, "#phi", "Centrality", "dN^{Jets}/d#phi");
214  AddHistogram2D<TH2D>(Form("hJetEta%s", appendix), "Jet angular distribution #eta", "LEGO2", 100, -2.5, 2.5, fNumberOfCentralityBins, 0, 100, "#eta","Centrality","dN^{Jets}/d#eta");
215  AddHistogram2D<TH2D>(Form("hJetPhiPt%s", appendix), "Jet angular distribution #phi vs. p_{T}", "LEGO2", 180, 0., 2*TMath::Pi(), 400, -100., 300., "#phi", "p_{T, jet} (GeV/c)", "dN^{Jets}/d#phi dp_{T}");
216  AddHistogram2D<TH2D>(Form("hJetEtaPt%s", appendix), "Jet angular distribution #eta vs. p_{T}", "LEGO2", 100, -2.5, 2.5, 400, -100., 300., "#eta","p_{T, jet} (GeV/c)","dN^{Jets}/d#eta dp_{T}");
217  AddHistogram2D<TH2D>(Form("hJetPhiEta%s", appendix), "Jet angular distribution #phi/#eta", "COLZ", 180, 0., 2*TMath::Pi(), 100, -2.5, 2.5, "#phi", "#eta", "dN^{Jets}/d#phi d#eta");
218  AddHistogram2D<TH2D>(Form("hJetArea%s", appendix), "Jet area", "LEGO2", 200, 0., 2., fNumberOfCentralityBins, 0, 100, "Jet A", "Centrality", "dN^{Jets}/dA");
219  AddHistogram2D<TH2D>(Form("hJetAreaPt%s", appendix), "Jet area vs. p_{T}", "LEGO2", 200, 0., 2., 400, -100., 300., "Jet A", "p_{T, jet} (GeV/c)", "dN^{Jets}/dA dp_{T}");
220  AddHistogram2D<TH2D>(Form("hJetPtLeadingHadron%s", appendix), "Jet leading hadron p_{T} distribution vs. jet p_{T}", "", 300, 0., 300., 300, 0., 300., "p_{T, jet} (GeV/c)", "p_{T,lead had} (GeV/c)", "dN^{Jets}/dp_{T}dp_{T,had}");
221 
222  // Leading/subleading ...
223  AddHistogram2D<TH2D>(Form("hLeadingJetPtRaw%s", appendix), "Jets p_{T} distribution (no bgrd. corr.)", "", 300, 0., 300., fNumberOfCentralityBins, 0, 100, "p_{T, jet} (GeV/c)", "Centrality", "dN^{Jets}/dp_{T}");
224  AddHistogram2D<TH2D>(Form("hLeadingJetPt%s", appendix), "Jets p_{T} distribution (background subtracted)", "", 400, -100., 300., fNumberOfCentralityBins, 0, 100, "p_{T, jet} (GeV/c)", "Centrality", "dN^{Jets}/dp_{T}");
225  AddHistogram2D<TH2D>(Form("hLeadingJetPhi%s", appendix), "Jet angular distribution #phi", "LEGO2", 180, 0., 2*TMath::Pi(), fNumberOfCentralityBins, 0, 100, "#phi", "Centrality", "dN^{Jets}/d#phi");
226  AddHistogram2D<TH2D>(Form("hLeadingJetEta%s", appendix), "Jet angular distribution #eta", "LEGO2", 100, -2.5, 2.5, fNumberOfCentralityBins, 0, 100, "#eta","Centrality","dN^{Jets}/d#eta");
227  AddHistogram2D<TH2D>(Form("hLeadingJetPhiPt%s", appendix), "Jet angular distribution #phi vs. p_{T}", "LEGO2", 180, 0., 2*TMath::Pi(), 400, -100., 300., "#phi", "p_{T, jet} (GeV/c)", "dN^{Jets}/d#phi dp_{T}");
228  AddHistogram2D<TH2D>(Form("hLeadingJetEtaPt%s", appendix), "Jet angular distribution #eta vs. p_{T}", "LEGO2", 100, -2.5, 2.5, 400, -100., 300., "#eta","p_{T, jet} (GeV/c)","dN^{Jets}/d#eta dp_{T}");
229  AddHistogram2D<TH2D>(Form("hLeadingJetPhiEta%s", appendix), "Jet angular distribution #phi/#eta", "COLZ", 180, 0., 2*TMath::Pi(), 100, -2.5, 2.5, "#phi", "#eta", "dN^{Jets}/d#phi d#eta");
230  AddHistogram2D<TH2D>(Form("hLeadingJetArea%s", appendix), "Jet area", "LEGO2", 200, 0., 2., fNumberOfCentralityBins, 0, 100, "Jet A", "Centrality", "dN^{Jets}/dA");
231  AddHistogram2D<TH2D>(Form("hLeadingJetAreaPt%s", appendix), "Jet area vs. p_{T}", "LEGO2", 200, 0., 2., 400, -100., 300., "Jet A", "p_{T, jet} (GeV/c)", "dN^{Jets}/dA dp_{T}");
232  AddHistogram2D<TH2D>(Form("hLeadingJetPtLeadingHadron%s", appendix), "Jet leading hadron p_{T} distribution vs. jet p_{T}", "", 300, 0., 300., 300, 0., 300., "p_{T, jet} (GeV/c)", "p_{T,lead had} (GeV/c)", "dN^{Jets}/dp_{T}dp_{T,had}");
233 
234  AddHistogram2D<TH2D>(Form("hSubleadingJetPtRaw%s", appendix), "Jets p_{T} distribution (no bgrd. corr.)", "", 300, 0., 300., fNumberOfCentralityBins, 0, 100, "p_{T, jet} (GeV/c)", "Centrality", "dN^{Jets}/dp_{T}");
235  AddHistogram2D<TH2D>(Form("hSubleadingJetPt%s", appendix), "Jets p_{T} distribution (background subtracted)", "", 400, -100., 300., fNumberOfCentralityBins, 0, 100, "p_{T, jet} (GeV/c)", "Centrality", "dN^{Jets}/dp_{T}");
236  AddHistogram2D<TH2D>(Form("hSubleadingJetPhi%s", appendix), "Jet angular distribution #phi", "LEGO2", 180, 0., 2*TMath::Pi(), fNumberOfCentralityBins, 0, 100, "#phi", "Centrality", "dN^{Jets}/d#phi");
237  AddHistogram2D<TH2D>(Form("hSubleadingJetEta%s", appendix), "Jet angular distribution #eta", "LEGO2", 100, -2.5, 2.5, fNumberOfCentralityBins, 0, 100, "#eta","Centrality","dN^{Jets}/d#eta");
238  AddHistogram2D<TH2D>(Form("hSubleadingJetPhiPt%s", appendix), "Jet angular distribution #phi vs. p_{T}", "LEGO2", 180, 0., 2*TMath::Pi(), 400, -100., 300., "#phi", "p_{T, jet} (GeV/c)", "dN^{Jets}/d#phi dp_{T}");
239  AddHistogram2D<TH2D>(Form("hSubleadingJetEtaPt%s", appendix), "Jet angular distribution #eta vs. p_{T}", "LEGO2", 100, -2.5, 2.5, 400, -100., 300., "#eta","p_{T, jet} (GeV/c)","dN^{Jets}/d#eta dp_{T}");
240  AddHistogram2D<TH2D>(Form("hSubleadingJetPhiEta%s", appendix), "Jet angular distribution #phi/#eta", "COLZ", 180, 0., 2*TMath::Pi(), 100, -2.5, 2.5, "#phi", "#eta", "dN^{Jets}/d#phi d#eta");
241  AddHistogram2D<TH2D>(Form("hSubleadingJetArea%s", appendix), "Jet area", "LEGO2", 200, 0., 2., fNumberOfCentralityBins, 0, 100, "Jet A", "Centrality", "dN^{Jets}/dA");
242  AddHistogram2D<TH2D>(Form("hSubleadingJetAreaPt%s", appendix), "Jet area vs. p_{T}", "LEGO2", 200, 0., 2., 400, -100., 300., "Jet A", "p_{T, jet} (GeV/c)", "dN^{Jets}/dA dp_{T}");
243  AddHistogram2D<TH2D>(Form("hSubleadingJetPtLeadingHadron%s", appendix), "Jet leading hadron p_{T} distribution vs. jet p_{T}", "", 300, 0., 300., 300, 0., 300., "p_{T, jet} (GeV/c)", "p_{T,lead had} (GeV/c)", "dN^{Jets}/dp_{T}dp_{T,had}");
244 
245  AddHistogram2D<TH2D>(Form("hJetConstituentPt_Cent0_100%s", appendix), "Jet constituent p_{T} distribution vs. jet p_T (background subtracted)", "", 400, -100., 300., 300, 0., 300., "p_{T, jet} (GeV/c)", "p_{T, track} (GeV/c)", "dN^{Tracks}/d^{2}p_{T}");
246  AddHistogram2D<TH2D>(Form("hJetConstituentPt_Cent0_10%s", appendix), "Jet constituent p_{T} distribution vs. jet p_T (background subtracted), 0-10 centrality", "", 400, -100., 300., 300, 0., 300., "p_{T, jet} (GeV/c)", "p_{T, track} (GeV/c)", "dN^{Tracks}/d^{2}p_{T}");
247  AddHistogram2D<TH2D>(Form("hJetConstituentPt_Cent10_30%s", appendix), "Jet constituent p_{T} distribution vs. jet p_T (background subtracted), 10-30 centrality", "", 400, -100., 300., 300, 0., 300., "p_{T, jet} (GeV/c)", "p_{T, track} (GeV/c)", "dN^{Tracks}/d^{2}p_{T}");
248  AddHistogram2D<TH2D>(Form("hJetConstituentPt_Cent30_50%s", appendix), "Jet constituent p_{T} distribution vs. jet p_T (background subtracted), 30-50 centrality", "", 400, -100., 300., 300, 0., 300., "p_{T, jet} (GeV/c)", "p_{T, track} (GeV/c)", "dN^{Tracks}/d^{2}p_{T}");
249  AddHistogram2D<TH2D>(Form("hJetConstituentPt_Cent50_90%s", appendix), "Jet constituent p_{T} distribution vs. jet p_T (background subtracted), 50-90 centrality", "", 400, -100., 300., 300, 0., 300., "p_{T, jet} (GeV/c)", "p_{T, track} (GeV/c)", "dN^{Tracks}/d^{2}p_{T}");
250 
251  AddHistogram2D<TH2D>(Form("hJetConstituentPt_Cent0_100_FilterBit%i%s", fConstPtFilterBit, appendix), "Jet constituent p_{T} distribution vs. jet p_T (background subtracted)", "", 400, -100., 300., 300, 0., 300., "p_{T, jet} (GeV/c)", "p_{T, track} (GeV/c)", "dN^{Tracks}/d^{2}p_{T}");
252  AddHistogram2D<TH2D>(Form("hJetConstituentPt_Cent0_10_FilterBit%i%s", fConstPtFilterBit, appendix), "Jet constituent p_{T} distribution vs. jet p_T (background subtracted), 0-10 centrality", "", 400, -100., 300., 300, 0., 300., "p_{T, jet} (GeV/c)", "p_{T, track} (GeV/c)", "dN^{Tracks}/d^{2}p_{T}");
253  AddHistogram2D<TH2D>(Form("hJetConstituentPt_Cent10_30_FilterBit%i%s", fConstPtFilterBit, appendix), "Jet constituent p_{T} distribution vs. jet p_T (background subtracted), 10-30 centrality", "", 400, -100., 300., 300, 0., 300., "p_{T, jet} (GeV/c)", "p_{T, track} (GeV/c)", "dN^{Tracks}/d^{2}p_{T}");
254  AddHistogram2D<TH2D>(Form("hJetConstituentPt_Cent30_50_FilterBit%i%s", fConstPtFilterBit, appendix), "Jet constituent p_{T} distribution vs. jet p_T (background subtracted), 30-50 centrality", "", 400, -100., 300., 300, 0., 300., "p_{T, jet} (GeV/c)", "p_{T, track} (GeV/c)", "dN^{Tracks}/d^{2}p_{T}");
255  AddHistogram2D<TH2D>(Form("hJetConstituentPt_Cent50_90_FilterBit%i%s", fConstPtFilterBit, appendix), "Jet constituent p_{T} distribution vs. jet p_T (background subtracted), 50-90 centrality", "", 400, -100., 300., 300, 0., 300., "p_{T, jet} (GeV/c)", "p_{T, track} (GeV/c)", "dN^{Tracks}/d^{2}p_{T}");
256 
257  AddHistogram2D<TH2D>(Form("hJetConstituentCount_Cent0_100%s", appendix), "Jet constituent count vs. jet p_T (background subtracted)", "", 400, -100., 300., 200, 0., 200., "p_{T, jet} (GeV/c)", "Count", "dN^{Jets}/dNdp_{T}");
258  AddHistogram2D<TH2D>(Form("hJetConstituentCount_Cent0_10%s", appendix), "Jet constituent count vs. jet p_T (background subtracted), 0-10 centrality", "", 400, -100., 300., 200, 0., 200., "p_{T, jet} (GeV/c)", "Count", "dN^{Jets}/dNdp_{T}");
259  }
260 
261  // Embedding plots
262  if(fJetOutputMode == 4 || fJetOutputMode == 5)
263  {
264  Double_t maxRatio = 1.;
266  maxRatio = 2.;
267 
268  for(Int_t i = -1; i<static_cast<Int_t>(fJetEmbeddingCuts.size()); i++)
269  {
270  const char* appendix = "";
271  if(i>-1)
272  {
274  appendix = Form("_%s", currentCut.fCutName.Data());
275 
276  // Don't double-add cuts
277  if( static_cast<TH1*>(fOutput->FindObject(Form("hEmbeddingDeltaR%s", appendix))) )
278  continue;
279  }
280  AddHistogram2D<TH2D>(Form("hEmbeddingDeltaR%s", appendix), "Matched jet #Delta R distribution", "", 200, -50., 150., 100, 0, 1.0, "p_{T, jet} (GeV/c)", "#Delta R", "dN^{Matched}/dp_{T}dR");
281  AddHistogram2D<TH2D>(Form("hEmbeddingDeltaEta%s", appendix), "Matched jet #Delta #eta distribution", "", 200, -50., 150., 100, -1.0, 1.0, "p_{T, jet} (GeV/c)", "#Delta #eta", "dN^{Matched}/dp_{T}d#eta");
282  AddHistogram2D<TH2D>(Form("hEmbeddingDeltaPhi%s", appendix), "Matched jet #Delta #phi distribution", "", 200, -50., 150., 100, -1.0, 1.0, "p_{T, jet} (GeV/c)", "#Delta #phi", "dN^{Matched}/dp_{T}d#phi");
283  AddHistogram2D<TH2D>(Form("hEmbeddingDeltaPt%s", appendix), "Matched jet #Delta p_{T} distribution", "", 200, -50., 150., 300, -150.0, 150.0, "p_{T, jet} (GeV/c)", "#Delta p_{T}", "dN^{Matched}/dp_{T}dp_{T}");
284  AddHistogram1D<TH1D>(Form("hEmbeddingJetPt%s", appendix), "Embedded jets p_{T} distribution", "", 200, -50., 150., "p_{T, jet} (GeV/c)", "dN/dp_{T}");
285  AddHistogram2D<TH2D>(Form("hEmbeddingJetPhiEta%s", appendix), "Embedded jet angular distribution #phi/#eta", "COLZ", 180, 0., 2*TMath::Pi(), 100, -2.5, 2.5, "#phi", "#eta", "dN^{Jets}/d#phi d#eta");
286 
288  {
289  AddHistogram3D<TH3D>(Form("hEmbeddingPtCorr010%s", appendix), "Matched jet p_{T} distributions (0-10% centrality)", "", 150, 0., 150., 150, 0., 150., 100, 0., maxRatio, "p_{T, MC jet} (GeV/c)", "p_{T, emb} (GeV/c)", "% MC");
290  AddHistogram3D<TH3D>(Form("hEmbeddingPtCorr1030%s", appendix), "Matched jet p_{T} distributions (10-30% centrality)", "", 150, 0., 150., 150, 0., 150., 100, 0., maxRatio, "p_{T, MC jet} (GeV/c)", "p_{T, emb} (GeV/c)", "% MC");
291  AddHistogram3D<TH3D>(Form("hEmbeddingPtCorr3050%s", appendix), "Matched jet p_{T} distributions (30-50% centrality)", "", 150, 0., 150., 150, 0., 150., 100, 0., maxRatio, "p_{T, MC jet} (GeV/c)", "p_{T, emb} (GeV/c)", "% MC");
292  AddHistogram3D<TH3D>(Form("hEmbeddingPtCorr5090%s", appendix), "Matched jet p_{T} distributions (50-90% centrality)", "", 150, 0., 150., 150, 0., 150., 100, 0., maxRatio, "p_{T, MC jet} (GeV/c)", "p_{T, emb} (GeV/c)", "% MC");
293  }
294  else
295  {
296  AddHistogram2D<TH2D>(Form("hEmbeddingPtCorr010_Above20%s", appendix), "Matched jet p_{T} distributions, MC ratio > 20% (0-10% centrality)", "", 150, 0., 150., 150, 0., 150.,"p_{T, MC jet} (GeV/c)", "p_{T, emb} (GeV/c)", "% MC");
297  AddHistogram2D<TH2D>(Form("hEmbeddingPtCorr1030_Above20%s", appendix),"Matched jet p_{T} distributions, MC ratio > 20% (10-30% centrality)","", 150, 0., 150., 150, 0., 150.,"p_{T, MC jet} (GeV/c)", "p_{T, emb} (GeV/c)", "% MC");
298  AddHistogram2D<TH2D>(Form("hEmbeddingPtCorr3050_Above20%s", appendix),"Matched jet p_{T} distributions, MC ratio > 20% (30-50% centrality)","", 150, 0., 150., 150, 0., 150.,"p_{T, MC jet} (GeV/c)", "p_{T, emb} (GeV/c)", "% MC");
299  AddHistogram2D<TH2D>(Form("hEmbeddingPtCorr5090_Above20%s", appendix),"Matched jet p_{T} distributions, MC ratio > 20% (50-90% centrality)","", 150, 0., 150., 150, 0., 150.,"p_{T, MC jet} (GeV/c)", "p_{T, emb} (GeV/c)", "% MC");
300  }
301  }
302 
304  {
305  AddHistogram3D<TH3D>("hEmbeddingPtCorr010", "Matched jet p_{T} distributions (0-10% centrality)", "", 150, 0., 150., 150, 0., 150., 100, 0., maxRatio, "p_{T, MC jet} (GeV/c)", "p_{T, emb} (GeV/c)", "% MC");
306  AddHistogram3D<TH3D>("hEmbeddingPtCorr1030", "Matched jet p_{T} distributions (10-30% centrality)", "", 150, 0., 150., 150, 0., 150., 100, 0., maxRatio, "p_{T, MC jet} (GeV/c)", "p_{T, emb} (GeV/c)", "% MC");
307  AddHistogram3D<TH3D>("hEmbeddingPtCorr3050", "Matched jet p_{T} distributions (30-50% centrality)", "", 150, 0., 150., 150, 0., 150., 100, 0., maxRatio, "p_{T, MC jet} (GeV/c)", "p_{T, emb} (GeV/c)", "% MC");
308  AddHistogram3D<TH3D>("hEmbeddingPtCorr5090", "Matched jet p_{T} distributions (50-90% centrality)", "", 150, 0., 150., 150, 0., 150., 100, 0., maxRatio, "p_{T, MC jet} (GeV/c)", "p_{T, emb} (GeV/c)", "% MC");
309  }
310  }
311 
312  // Random cone plots, background, ...
313  AddHistogram2D<TH2D>("hRandomConePt", "Random cone p_{T} distribution", "", 400, -100., 300., fNumberOfCentralityBins, 0, 100, "p_{T, cone} (GeV/c)", "Centrality", "dN^{Tracks}/dp_{T}");
314  AddHistogram2D<TH2D>("hRandomConePtCut3GeV", "Random cone p_{T} distribution, cut p_{T} > 3 GeV/c", "", 400, -100., 300., fNumberOfCentralityBins, 0, 100, "p_{T, cone} (GeV/c)", "Centrality", "dN^{Tracks}/dp_{T}");
315  AddHistogram2D<TH2D>("hRandomConeRawPt", "Random cone p_{T} distribution (no bgrd. correction)", "", 300, 0., 300., fNumberOfCentralityBins, 0, 100, "p_{T, cone} (GeV/c)", "Centrality", "dN^{Tracks}/dp_{T}");
316  AddHistogram2D<TH2D>("hRandomConeRawPtCut3GeV", "Random cone p_{T} distribution (no bgrd. correction), cut p_{T} > 3 GeV/c", "", 300, 0., 300., fNumberOfCentralityBins, 0, 100, "p_{T, cone} (GeV/c)", "Centrality", "dN^{Tracks}/dp_{T}");
317 
318  AddHistogram2D<TH2D>("hTrackCount", "Number of tracks in acceptance vs. centrality", "LEGO2", 500, 0., 5000., fNumberOfCentralityBins, 0, 100, "N tracks","Centrality", "dN^{Events}/dN^{Tracks}");
319  AddHistogram2D<TH2D>("hJetCount", "Number of jets in acceptance vs. centrality", "LEGO2", 100, 0., 100., fNumberOfCentralityBins, 0, 100, "N Jets","Centrality", "dN^{Events}/dN^{Jets}");
320  AddHistogram2D<TH2D>("hBackgroundPt", "Background p_{T} distribution", "", 150, 0., 150., fNumberOfCentralityBins, 0, 100, "Background p_{T} (GeV/c)", "Centrality", "dN^{Events}/dp_{T}");
321 
322 
323  for(Int_t i = -1; i<static_cast<Int_t>(fJetEmbeddingCuts.size()); i++)
324  {
325  const char* appendix = "";
326  if(i>-1)
327  {
329  appendix = Form("_%s", currentCut.fCutName.Data());
330 
331  // Don't double-add cuts
332  if( static_cast<TH1*>(fOutput->FindObject(Form("hBackgroundPtJetPt_Cent0_100%s", appendix))) )
333  continue;
334  }
335  AddHistogram2D<TH2D>(Form("hBackgroundPtJetPt_Cent0_100%s", appendix), "Background p_{T} distribution vs. jet p_{T}", "", 150, 0., 150., 400, -100., 300., "Background p_{T} (GeV/c)", "Jet p_{T} (GeV/c)", "dN^{Events}/dp_{T}");
336  AddHistogram2D<TH2D>(Form("hBackgroundPtJetPt_Cent0_10%s", appendix), "Background p_{T} distribution vs. jet p_{T}", "", 150, 0., 150., 400, -100., 300., "Background p_{T} (GeV/c)", "Jet p_{T} (GeV/c)", "dN^{Events}/dp_{T}");
337  AddHistogram2D<TH2D>(Form("hBackgroundPtConstCount_Cent0_100%s", appendix), "Background p_{T} distribution vs. const. count", "", 150, 0., 150., 200, 0., 200., "Background p_{T} (GeV/c)", "Count", "dN^{Events}/dp_{T}");
338  AddHistogram2D<TH2D>(Form("hBackgroundPtConstCount_Cent0_10%s", appendix), "Background p_{T} distribution vs. const. count", "", 150, 0., 150., 200, 0., 200., "Background p_{T} (GeV/c)", "Count", "dN^{Events}/dp_{T}");
339  }
340 
341  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
342 }
343 
344 
345 //________________________________________________________________________
347 
349 
350 
351  // ### Add the jets as basic correlation particles to the event
352  // This output object carries all accepted jets
353  if (!(fInputEvent->FindListObject(Form("%s", fJetParticleArrayName.Data()))))
354  {
355  fJetsOutput.push_back(new TClonesArray("AliPicoTrack"));
356  fJetsOutput.at(0)->SetName(fJetParticleArrayName.Data());
357  fInputEvent->AddObject(fJetsOutput.at(0));
358  }
359  else
360  AliFatal(Form("%s: Object with name %s already in event!", GetName(), Form("%s", fJetParticleArrayName.Data())));
361 
362  // These output objects carry all jets that pass certain cuts
363  if( (fJetOutputMode==4 || fJetOutputMode==5) )
364  {
365  // before, check if all given output names are OK
366  for(Int_t i = 0; i<fJetEmbeddingCuts.size(); i++)
367  if (fInputEvent->FindListObject(Form("%s", fJetEmbeddingCuts.at(i).fOutputName.Data())))
368  AliFatal(Form("%s: Object with name %s already in event!", GetName(), Form("%s", fJetEmbeddingCuts.at(i).fOutputName.Data())));
369 
370  for(Int_t i = 0; i<fJetEmbeddingCuts.size(); i++)
371  {
372  // If the cut demands a new output stream, add it
373  if (!fInputEvent->FindListObject(Form("%s", fJetEmbeddingCuts.at(i).fOutputName.Data())))
374  {
375  fJetsOutput.push_back(new TClonesArray("AliPicoTrack"));
376  fJetsOutput.at(fJetsOutput.size()-1)->SetName(fJetEmbeddingCuts.at(i).fOutputName.Data());
377  fInputEvent->AddObject(fJetsOutput.at(fJetsOutput.size()-1));
378  fJetEmbeddingCuts.at(i).fArrayIndex = fJetsOutput.size()-1;
379 
380  // Set the array indices for all cuts that use this output stream
381  for(Int_t j = 0; j<fJetEmbeddingCuts.size(); j++)
382  {
383  if(fJetEmbeddingCuts.at(j).fArrayIndex != -1)
384  continue;
385  if(fJetEmbeddingCuts.at(j).fOutputName != fJetEmbeddingCuts.at(i).fOutputName)
386  continue;
387  fJetEmbeddingCuts.at(j).fArrayIndex = fJetEmbeddingCuts.at(i).fArrayIndex;
388  }
389  }
390  }
391  }
392 
393  // ##############################################
394  // ##############################################
395 
396  // ### Prepare the track tree
398  {
399  fTracksTree = new TTree("ExtractedTracks", "ExtractedTracks");
400  fTracksTree->Branch("Kinematics", "AliBasicParticle", &fTreeBufferTrack, 1000);
401  fTracksTree->Branch("PID", "AliAODPid", &fTreeBufferPID, 1000);
402  fTracksTree->Branch("PDG",&fTreeBufferPDG,"a/I");
403  fOutput->Add(fTracksTree);
404  }
405 
406  // ### Add the tracks as basic correlation particles to the event (optional)
407  if(fTrackParticleArrayName != "")
408  {
409  if (!(fInputEvent->FindListObject(Form("%s", fTrackParticleArrayName.Data()))))
410  {
411  fTracksOutput = new TClonesArray("AliPicoTrack");
412  fTracksOutput->SetName(fTrackParticleArrayName.Data());
413  fInputEvent->AddObject(fTracksOutput);
414  }
415  else
416  AliFatal(Form("%s: Object with name %s already in event!", GetName(), Form("%s", fTrackParticleArrayName.Data())));
417  }
418 
419  // ### Import jets for embedding (optional)
420  if(fJetEmbeddingArrayName != "")
421  {
422  fJetEmbeddingArray = static_cast<TClonesArray*>(InputEvent()->FindListObject(Form("%s", fJetEmbeddingArrayName.Data())));
423  if(!fJetEmbeddingArray)
424  AliFatal(Form("Importing jets for embedding failed! Array '%s' not found!", fJetEmbeddingArrayName.Data()));
425  }
426  else if(fJetOutputMode==4 || fJetOutputMode==5)
427  AliFatal(Form("fJetEmbeddingArrayName must be set in jet output mode 4 or 5."));
428 
429  // ### Import veto jets for matching (optional)
430  if(fJetVetoArrayName != "")
431  {
432  fJetVetoArray = static_cast<TClonesArray*>(InputEvent()->FindListObject(Form("%s", fJetVetoArrayName.Data())));
433  if(!fJetVetoArray)
434  AliFatal(Form("Importing jets for veto failed! Array '%s' not found!", fJetVetoArrayName.Data()));
435  }
436 }
437 
438 //________________________________________________________________________
440 {
441  if(fJetOutputMode==3) // output leading&subleading jet
442  {
443  if((jet!=fLeadingJet) && (jet!=fSubleadingJet))
444  return kFALSE;
445  }
446  else if(fJetOutputMode==1) // output the leading jet
447  {
448  if(jet!=fLeadingJet)
449  return kFALSE;
450  }
451  else if(fJetOutputMode==2) // output the subleading jet
452  {
453  if(jet!=fSubleadingJet)
454  return kFALSE;
455  }
456  else if(fJetOutputMode==6)
457  {
459  return kFALSE;
460  }
461 
462  if(fJetOutputMode==4) // matching jets only
463  return (std::find(fMatchedJets.begin(), fMatchedJets.end(), jet) != fMatchedJets.end());
464  else if(fJetOutputMode==5) // non-matching jets only
465  return (std::find(fMatchedJets.begin(), fMatchedJets.end(), jet) == fMatchedJets.end());
466 
467  return kTRUE;
468 }
469 
470 //________________________________________________________________________
472 {
473  TString appendix("");
474  if(cutName)
475  appendix = Form("_%s", cutName);
476 
477  // All jets
478  FillHistogram(Form("hJetPtRaw%s", appendix.Data()), jet->Pt(), fCent);
479  FillHistogram(Form("hJetPt%s", appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), fCent);
480  FillHistogram(Form("hJetPhi%s", appendix.Data()), jet->Phi(), fCent);
481  FillHistogram(Form("hJetEta%s", appendix.Data()), jet->Eta(), fCent);
482  FillHistogram(Form("hJetEtaPt%s", appendix.Data()), jet->Eta(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
483  FillHistogram(Form("hJetPhiPt%s", appendix.Data()), jet->Phi(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
484  FillHistogram(Form("hJetPhiEta%s", appendix.Data()), jet->Phi(), jet->Eta());
485  FillHistogram(Form("hJetArea%s", appendix.Data()), jet->Area(), fCent);
486  FillHistogram(Form("hJetAreaPt%s", appendix.Data()), jet->Area(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
487  FillHistogram(Form("hJetPtLeadingHadron%s", appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), fJetsCont->GetLeadingHadronPt(jet));
488 
489  FillHistogram(Form("hBackgroundPtJetPt_Cent0_100%s", appendix.Data()), fJetsCont->GetRhoVal(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
490  if( (fCent >= 0) && (fCent < 10) )
491  FillHistogram(Form("hBackgroundPtJetPt_Cent0_10%s", appendix.Data()), fJetsCont->GetRhoVal(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
492 
493  // Leading jet plots
494  if(jet==fLeadingJet)
495  {
496  FillHistogram(Form("hLeadingJetPtRaw%s", appendix.Data()), jet->Pt(), fCent);
497  FillHistogram(Form("hLeadingJetPt%s", appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), fCent);
498  FillHistogram(Form("hLeadingJetPhi%s", appendix.Data()), jet->Phi(), fCent);
499  FillHistogram(Form("hLeadingJetEta%s", appendix.Data()), jet->Eta(), fCent);
500  FillHistogram(Form("hLeadingJetEtaPt%s", appendix.Data()), jet->Eta(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
501  FillHistogram(Form("hLeadingJetPhiPt%s", appendix.Data()), jet->Phi(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
502  FillHistogram(Form("hLeadingJetPhiEta%s", appendix.Data()), jet->Phi(), jet->Eta());
503  FillHistogram(Form("hLeadingJetArea%s", appendix.Data()), jet->Area(), fCent);
504  FillHistogram(Form("hLeadingJetAreaPt%s", appendix.Data()), jet->Area(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
505  FillHistogram(Form("hLeadingJetPtLeadingHadron%s", appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), fJetsCont->GetLeadingHadronPt(jet));
506  }
507 
508  // Subleading jet plot
509  else if(jet==fSubleadingJet)
510  {
511  FillHistogram(Form("hSubleadingJetPtRaw%s", appendix.Data()), jet->Pt(), fCent);
512  FillHistogram(Form("hSubleadingJetPt%s", appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), fCent);
513  FillHistogram(Form("hSubleadingJetPhi%s", appendix.Data()), jet->Phi(), fCent);
514  FillHistogram(Form("hSubleadingJetEta%s", appendix.Data()), jet->Eta(), fCent);
515  FillHistogram(Form("hSubleadingJetEtaPt%s", appendix.Data()), jet->Eta(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
516  FillHistogram(Form("hSubleadingJetPhiPt%s", appendix.Data()), jet->Phi(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
517  FillHistogram(Form("hSubleadingJetPhiEta%s", appendix.Data()), jet->Phi(), jet->Eta());
518  FillHistogram(Form("hSubleadingJetArea%s", appendix.Data()), jet->Area(), fCent);
519  FillHistogram(Form("hSubleadingJetAreaPt%s", appendix.Data()), jet->Area(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
520  FillHistogram(Form("hSubleadingJetPtLeadingHadron%s", appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), fJetsCont->GetLeadingHadronPt(jet));
521  }
522 
523  // ####### Jet constituent plots
524  Int_t nProcessedTracks = 0;
525  for(Int_t i = 0; i < jet->GetNumberOfTracks(); i++)
526  {
527  AliVParticle* constituent = static_cast<AliVParticle*>(jet->TrackAt(i, fTracksCont->GetArray()));
528  if(!constituent)
529  continue;
530 
531  // Check whether to discard this track
532  if((!fUseMCConstituents) && (constituent->GetLabel() >= 10000)) // is from MC
533  continue;
534  if((!fUseDataConstituents) && (constituent->GetLabel() < 10000)) // is from data
535  continue;
536 
537  nProcessedTracks++;
538 
539  Bool_t filterConditionFulfilled = kFALSE;
540  AliAODTrack* aodTrack = static_cast<AliAODTrack*>(constituent);
541  if (aodTrack)
542  filterConditionFulfilled = aodTrack->TestFilterBit(fConstPtFilterBit);
543 
544 
545  // Fill jet constituent plots
546  FillHistogram(Form("hJetConstituentPt_Cent0_100%s", appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), constituent->Pt());
547  if(filterConditionFulfilled)
548  FillHistogram(Form("hJetConstituentPt_Cent0_100_FilterBit%i%s", fConstPtFilterBit, appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), constituent->Pt());
549  if( (fCent >= 0) && (fCent < 10) )
550  {
551  FillHistogram(Form("hJetConstituentPt_Cent0_10%s", appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), constituent->Pt());
552  if(filterConditionFulfilled)
553  FillHistogram(Form("hJetConstituentPt_Cent0_10_FilterBit%i%s", fConstPtFilterBit, appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), constituent->Pt());
554  }
555  else if( (fCent >= 10) && (fCent < 30) )
556  {
557  FillHistogram(Form("hJetConstituentPt_Cent10_30%s", appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), constituent->Pt());
558  if(filterConditionFulfilled)
559  FillHistogram(Form("hJetConstituentPt_Cent10_30_FilterBit%i%s", fConstPtFilterBit, appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), constituent->Pt());
560  }
561  else if( (fCent >= 30) && (fCent < 50) )
562  {
563  FillHistogram(Form("hJetConstituentPt_Cent30_50%s", appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), constituent->Pt());
564  if(filterConditionFulfilled)
565  FillHistogram(Form("hJetConstituentPt_Cent30_50_FilterBit%i%s", fConstPtFilterBit, appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), constituent->Pt());
566  }
567  else if( (fCent >= 50) && (fCent < 90) )
568  {
569  FillHistogram(Form("hJetConstituentPt_Cent50_90%s", appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), constituent->Pt());
570  if(filterConditionFulfilled)
571  FillHistogram(Form("hJetConstituentPt_Cent50_90_FilterBit%i%s", fConstPtFilterBit, appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), constituent->Pt());
572  }
573  }
574 
575  FillHistogram(Form("hJetConstituentCount_Cent0_100%s", appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), nProcessedTracks);
576  if( (fCent >= 0) && (fCent < 10) )
577  FillHistogram(Form("hJetConstituentCount_Cent0_10%s", appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), nProcessedTracks);
578 
579  FillHistogram(Form("hBackgroundPtConstCount_Cent0_100%s", appendix.Data()), fJetsCont->GetRhoVal(), nProcessedTracks);
580  if( (fCent >= 0) && (fCent < 10) )
581  FillHistogram(Form("hBackgroundPtConstCount_Cent0_10%s", appendix.Data()), fJetsCont->GetRhoVal(), nProcessedTracks);
582 
583  // ####### Embedding plots
584  if( (fJetOutputMode == 4) || (fJetOutputMode == 5))
585  {
586  AliEmcalJet* refJet = GetReferenceJet(jet);
587  Double_t deltaPt = jet->Pt() - fJetsCont->GetRhoVal()*jet->Area() - refJet->Pt();
588  Double_t deltaEta = (jet->Eta()-refJet->Eta());
589  Double_t deltaPhi = TMath::Min(TMath::Abs(jet->Phi()-refJet->Phi()),TMath::TwoPi() - TMath::Abs(jet->Phi()-refJet->Phi()));
590  if(jet->Phi() < refJet->Phi())
591  deltaPhi = -deltaPhi;
592 
593  Double_t deltaR = TMath::Sqrt(deltaEta*deltaEta + deltaPhi*deltaPhi);
594  FillHistogram(Form("hEmbeddingDeltaR%s", appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), deltaR);
595  FillHistogram(Form("hEmbeddingDeltaEta%s", appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), deltaPhi);
596  FillHistogram(Form("hEmbeddingDeltaPhi%s", appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), deltaEta);
597  FillHistogram(Form("hEmbeddingDeltaPt%s", appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), deltaPt);
598  FillHistogram(Form("hEmbeddingJetPt%s", appendix.Data()), refJet->Pt());
599  FillHistogram(Form("hEmbeddingJetPhiEta%s", appendix.Data()), refJet->Phi(), refJet->Eta());
600 
601  // Only create 3D plots for each cut on demand
602  Double_t trackRatio = 0.;
603  Double_t ptRatio = 0.;
604  GetTrackMCRatios(jet, refJet, trackRatio, ptRatio);
605 
606  if(fCent >= 0 && fCent < 10)
607  {
608  if((appendix == "") || fJetEmbeddingCreatePtPlotPerCut)
609  FillHistogram3D(Form("hEmbeddingPtCorr010%s", appendix.Data()), refJet->Pt(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), ptRatio);
610  if(ptRatio >= 0.2)
611  FillHistogram(Form("hEmbeddingPtCorr010_Above20%s", appendix.Data()), refJet->Pt(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
612  }
613  else if (fCent >= 10 && fCent < 30)
614  {
615  if((appendix == "") || fJetEmbeddingCreatePtPlotPerCut)
616  FillHistogram3D(Form("hEmbeddingPtCorr1030%s", appendix.Data()), refJet->Pt(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), ptRatio);
617  if(ptRatio >= 0.2)
618  FillHistogram(Form("hEmbeddingPtCorr1030_Above20%s", appendix.Data()), refJet->Pt(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
619  }
620  else if (fCent >= 30 && fCent < 50)
621  {
622  if((appendix == "") || fJetEmbeddingCreatePtPlotPerCut)
623  FillHistogram3D(Form("hEmbeddingPtCorr3050%s", appendix.Data()), refJet->Pt(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), ptRatio);
624  if(ptRatio >= 0.2)
625  FillHistogram(Form("hEmbeddingPtCorr3050_Above20%s", appendix.Data()), refJet->Pt(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
626  }
627  else if (fCent >= 50 && fCent < 90)
628  {
629  if((appendix == "") || fJetEmbeddingCreatePtPlotPerCut)
630  FillHistogram3D(Form("hEmbeddingPtCorr5090%s", appendix.Data()), refJet->Pt(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), ptRatio);
631  if(ptRatio >= 0.2)
632  FillHistogram(Form("hEmbeddingPtCorr5090_Above20%s", appendix.Data()), refJet->Pt(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
633  }
634  }
635 }
636 
637 //________________________________________________________________________
639 {
640  FillHistogram("hTrackPt", track->Pt(), fCent);
641  FillHistogram("hTrackPhi", track->Phi(), fCent);
642  FillHistogram("hTrackEta", track->Eta(), fCent);
643  FillHistogram("hTrackEtaPt", track->Eta(), track->Pt());
644  FillHistogram("hTrackPhiPt", track->Phi(), track->Pt());
645  FillHistogram("hTrackPhiEta", track->Phi(), track->Eta());
646 }
647 
648 //________________________________________________________________________
650 {
651  Double_t tmpPt = jet->Pt() - fJetsCont->GetRhoVal()*jet->Area();
652  new ((*fJetsOutput.at(arrayIndex))[jetsAlreadyInArray]) AliPicoTrack(tmpPt, jet->Eta(), jet->Phi(), jet->Charge(), 0, 0);
653  jetsAlreadyInArray++;
654 }
655 
656 //________________________________________________________________________
658 {
659  if(fTrackParticleArrayName != "")
660  {
661  new ((*fTracksOutput)[fAcceptedTracks]) AliPicoTrack(track->Pt(), track->Eta(), track->Phi(), track->Charge(), 0, 0); // only Pt,Eta,Phi are interesting for correlations;
662  fAcceptedTracks++;
663  }
664 }
665 
666 
667 //________________________________________________________________________
669 {
670  // Only allow when we have aod tracks
671  AliAODTrack* aodtrack = dynamic_cast<AliAODTrack*>(track);
672  if(!aodtrack)
673  return;
674 
675  // Discard tracks according to their pT (below 20 GeV)
676  if(track->Pt() < 20.)
677  if(fRandom->Rndm() >= TMath::Power((1/20. * track->Pt()), fTrackExtractionPercentagePower))
678  return;
679 
680  // Create basic particle from track and extract pid object
681  fTreeBufferPID = aodtrack->GetDetPid();
682 
683  Int_t truthPID = 0;
684 
685  // Get truth values if we are on MC
686  TClonesArray* fTruthParticleArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("mcparticles"));
687  if(fTruthParticleArray)
688  {
689  for(Int_t i=0; i<fTruthParticleArray->GetEntries();i++)
690  {
691  AliAODMCParticle* mcParticle = dynamic_cast<AliAODMCParticle*>(fTruthParticleArray->At(i));
692  if(!mcParticle) continue;
693 
694  if (mcParticle->GetLabel() == aodtrack->GetLabel())
695  {
696  truthPID = mcParticle->PdgCode();
697  break;
698  }
699  }
700  }
701 
702  fTreeBufferPDG = truthPID;
703 
704  AliBasicParticle basicParticle(aodtrack->Eta(), aodtrack->Phi(), aodtrack->Pt(), aodtrack->Charge());
705  fTreeBufferTrack = &basicParticle;
706 
707  fTracksTree->Fill();
708 }
709 
710 //________________________________________________________________________
712 {
713  // Check jet pT threshold
715  return;
716 
717  // Discard jets statistically
718  if(fRandom->Rndm() >= fEventExtractionPercentage)
719  return;
720 
721  static Int_t numSavedEvents = 0;
722  numSavedEvents++;
723 
724 
725  AddHistogram2D<TH2D>(Form("Event%i", numSavedEvents), "Event display", "COLZ", 180, 0., 2*TMath::Pi(), 100, -2.5, 2.5, "#phi", "#eta", "dN^{Tracks}/d#phi d#eta");
726  fTracksCont->ResetCurrentID();
727  while(AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle()))
728  FillHistogram(Form("Event%i", numSavedEvents), track->Phi(), track->Eta(), track->Pt());
729 
730 }
731 
732 //________________________________________________________________________
734 {
736 
737  // Jet veto:
738  // 1. jet-by-jet mode: veto is active if jets overlaps with a suitable jet
739  // 2. other mode: veto is active if suitable jet is in sample
740  AliEmcalJet* vetoJet = 0;
741  if(!fJetVetoJetByJet)
742  vetoJet = GetLeadingVetoJet();
743 
744  // ####### Jet loop
745  fAcceptedJets = 0;
746  for(Int_t i = 0; i<fJetEmbeddingCuts.size(); i++)
747  fJetEmbeddingCuts.at(i).fAcceptedJets = 0;
748  fJetsCont->ResetCurrentID();
749  while(AliEmcalJet *jet = fJetsCont->GetNextAcceptJet())
750  {
751  if(!IsJetSelected(jet))
752  continue;
753 
754  // Plots + output jets regardless of embedding cut
755  FillHistogramsJets(jet, 0);
757 
758  // Plots + output jets for each embedding cut
759  if(fJetVetoJetByJet)
760  vetoJet = GetVetoJet(jet);
761  for(Int_t i = 0; i<fJetEmbeddingCuts.size(); i++)
762  {
763  AliChargedJetsHadronCFCuts* currentCut = &fJetEmbeddingCuts.at(i);
764  AliEmcalJet* refJet = GetReferenceJet(jet);
765  Double_t trackRatio = 0.;
766  Double_t ptRatio = 0.;
767  GetTrackMCRatios(jet, refJet, trackRatio, ptRatio);
768 
769  Double_t vetoJetPt = 0.;
770  if(vetoJet)
771  vetoJetPt = vetoJet->Pt() - vetoJet->Area()*fJetsCont->GetRhoVal();
772 
773  if(!currentCut->IsCutFulfilled(jet->Pt(), refJet->Pt(), fCent, ptRatio, vetoJetPt))
774  continue;
775 
776  FillHistogramsJets(jet, currentCut->fCutName.Data());
777  AddJetToOutputArray(jet, currentCut->fArrayIndex, currentCut->fAcceptedJets);
778  }
779  }
780 
781  // ####### Particle loop
782  fAcceptedTracks = 0;
783  fTracksCont->ResetCurrentID();
784  Int_t trackcount = 0;
785  while(AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle()))
786  {
787  // Check whether to discard this track
788  if((!fUseMCConstituents) && (track->GetLabel() >= 10000)) // is from MC
789  continue;
790  if((!fUseDataConstituents) && (track->GetLabel() < 10000)) // is from data
791  continue;
792 
793  // Track plots
794  FillHistogramsTracks(track);
795  // Add track to output array
796  trackcount++;
797  AddTrackToOutputArray(track);
799  AddTrackToTree(track);
800  }
801 
802  // Add event to output tree
804  AddEventToTree();
805 
806  // ######### Random cone sampling
807  for(Int_t iCone=0; iCone<fNumRandomConesPerEvent; iCone++)
808  {
809  // Throw random cone
810  Double_t tmpRandConeEta = fJetsCont->GetJetEtaMin() + fRandom->Rndm()*TMath::Abs(fJetsCont->GetJetEtaMax()-fJetsCont->GetJetEtaMin());
811  Double_t tmpRandConePhi = fRandom->Rndm()*TMath::TwoPi();
812  Double_t tmpRandConePt = 0;
813  Double_t tmpRandConePt3GeV = 0;
814  // Fill pT that is in cone
815  fTracksCont->ResetCurrentID();
816  while(AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle()))
817  {
818  // Check whether to discard this track
819  if((!fUseMCConstituents) && (track->GetLabel() >= 10000)) // is from MC
820  continue;
821  if((!fUseDataConstituents) && (track->GetLabel() < 10000)) // is from data
822  continue;
823  if(IsTrackInCone(track, tmpRandConeEta, tmpRandConePhi, fJetsCont->GetJetRadius()))
824  {
825  tmpRandConePt += track->Pt();
826  if (track->Pt() > 3.0)
827  tmpRandConePt3GeV += track->Pt();
828  }
829  }
830  // Fill histograms
831  FillHistogram("hRandomConePt", tmpRandConePt - fJetsCont->GetRhoVal()*fJetsCont->GetJetRadius()*fJetsCont->GetJetRadius()*TMath::Pi(), fCent);
832  FillHistogram("hRandomConePtCut3GeV", tmpRandConePt3GeV - fJetsCont->GetRhoVal()*fJetsCont->GetJetRadius()*fJetsCont->GetJetRadius()*TMath::Pi(), fCent);
833  FillHistogram("hRandomConeRawPt", tmpRandConePt, fCent);
834  FillHistogram("hRandomConeRawPtCut3GeV", tmpRandConePt3GeV, fCent);
835  }
836 
837  // ####### Event properties
838  FillHistogram("hBackgroundPt", fJetsCont->GetRhoVal(), fCent);
839  FillHistogram("hJetCount", fAcceptedJets, fCent);
840  FillHistogram("hTrackCount", trackcount, fCent);
841  // NOTE: It is possible to use fTracksCont->GetLeadingParticle() since we do not apply additional track cuts
842  AliVTrack* leadTrack = static_cast<AliVTrack*>(fTracksCont->GetLeadingParticle());
843  if(leadTrack)
844  {
845  FillHistogram("hLeadingTrackPt", leadTrack->Pt(), fCent);
846  FillHistogram("hLeadingTrackPhi", leadTrack->Phi(), fCent);
847  FillHistogram("hLeadingTrackEta", leadTrack->Eta(), fCent);
848  FillHistogram("hLeadingTrackPhiEta", leadTrack->Phi(), leadTrack->Eta());
849  }
850 
851  return kTRUE;
852 }
853 
854 //########################################################################
855 // HELPERS
856 //########################################################################
857 
858 //________________________________________________________________________
860 {
861  if(!fPythiaInfo)
862  AliError("fPythiaInfo object not available. Is it activated with SetGeneratePythiaInfoObject()?");
863 
864  Double_t bestMatchDeltaR1 = 999.;
865  Double_t bestMatchDeltaR2 = 999.;
866 
867  fJetsCont->ResetCurrentID();
868  while(AliEmcalJet *jet = fJetsCont->GetNextAcceptJet())
869  {
870  // Check via geometrical matching if jet is connected to the initial collision
871  Double_t deltaEta1 = TMath::Abs(jet->Eta()-fPythiaInfo->GetPartonEta6());
872  Double_t deltaEta2 = TMath::Abs(jet->Eta()-fPythiaInfo->GetPartonEta7());
873  Double_t deltaPhi1 = TMath::Min(TMath::Abs(jet->Phi()-fPythiaInfo->GetPartonPhi6()),TMath::TwoPi() - TMath::Abs(jet->Phi()-fPythiaInfo->GetPartonPhi6()));
874  Double_t deltaPhi2 = TMath::Min(TMath::Abs(jet->Phi()-fPythiaInfo->GetPartonPhi7()),TMath::TwoPi() - TMath::Abs(jet->Phi()-fPythiaInfo->GetPartonPhi7()));
875 
876  Double_t deltaR1 = TMath::Sqrt(deltaEta1*deltaEta1 + deltaPhi1*deltaPhi1);
877  Double_t deltaR2 = TMath::Sqrt(deltaEta2*deltaEta2 + deltaPhi2*deltaPhi2);
878 
879  if(deltaR1 < bestMatchDeltaR1)
880  {
881  bestMatchDeltaR1 = deltaR1;
883  }
884  if(deltaR2 < bestMatchDeltaR2)
885  {
886  bestMatchDeltaR2 = deltaR2;
888  }
889  }
890 
891  if(bestMatchDeltaR1 > 0.3)
893  if(bestMatchDeltaR2 > 0.3)
895 }
896 
897 //________________________________________________________________________
899 {
900  fMatchedJets.clear();
901  fMatchedJetsReference.clear();
902 
903  AliEmcalJet* jetLeading = 0;
904  AliEmcalJet* jetSubLeading = 0;
905 
907  GetLeadingJetsInArray(fJetEmbeddingArray, "", jetLeading, jetSubLeading);
908 
909  // ############ Search for all matches
910  for(Int_t i=0; i<TMath::Min(fJetEmbeddingArray->GetEntries(), fJetEmbeddingNumMatchedJets); i++)
911  {
912  AliEmcalJet* probeJet = 0;
913  // Extract leading/subleading
915  {
916  if(i==0)
917  probeJet = jetLeading;
918  else if(i==1)
919  probeJet = jetSubLeading;
920  }
921  else // extract more than 2 jets, e.g. all
922  probeJet = static_cast<AliEmcalJet*>(fJetEmbeddingArray->At(i));
923 
924  if(!probeJet)
925  continue;
926 
927  if(probeJet->Pt() < 0.001) // do not use ghosts
928  continue;
929 
930  AliEmcalJet* matchedJet = 0;
931  AliEmcalJet* matchedJetReference = 0;
932  Double_t bestMatchDeltaR = 999.;
933  fJetsCont->ResetCurrentID();
934  // Loop over all embedded jets to find the best match
935  while(AliEmcalJet* embeddedJet = fJetsCont->GetNextAcceptJet())
936  {
937  Double_t deltaEta = (embeddedJet->Eta()-probeJet->Eta());
938  Double_t deltaPhi = TMath::Min(TMath::Abs(embeddedJet->Phi()-probeJet->Phi()),TMath::TwoPi() - TMath::Abs(embeddedJet->Phi()-probeJet->Phi()));
939  Double_t deltaR = TMath::Sqrt(deltaEta*deltaEta + deltaPhi*deltaPhi);
940 
941  // Cut jets too far away
942  if (deltaR > fJetEmbeddingMaxDistance)
943  continue;
944 
945  // Search for the best match
946  if(deltaR < bestMatchDeltaR)
947  {
948  bestMatchDeltaR = deltaR;
949  matchedJet = embeddedJet;
950  matchedJetReference = probeJet;
951  }
952  }
953  // Put matched jet to a list
954  if(matchedJet && matchedJetReference)
955  {
956  fMatchedJets.push_back(matchedJet);
957  fMatchedJetsReference.push_back(matchedJetReference);
958  }
959  }
960 }
961 
962 
963 //________________________________________________________________________
965 {
966  Int_t tracksFromMC = 0;
967  Int_t tracksTotal = 0;
968  Double_t ptFromMC = 0;
969  Double_t ptTotal = 0;
970 
971  if(fJetEmbeddingUsePerTrackMCPercentage) // Calculate MC track percentage from tracks of the matched jet
972  {
973  for(Int_t j = 0; j < jet->GetNumberOfTracks(); j++)
974  {
975  AliVParticle* constituent = static_cast<AliVParticle*>(jet->TrackAt(j, fTracksCont->GetArray()));
976  if(!constituent)
977  continue;
978 
979  TClonesArray* mcArray = static_cast<TClonesArray*>(InputEvent()->FindListObject(Form("%s", fJetEmbeddingTrackArrayName.Data())));
980  Bool_t foundInMC = kFALSE;
981  for(Int_t k=0; k<mcJet->GetNumberOfTracks(); k++)
982  {
983  AliVParticle* mcConstituent = static_cast<AliVParticle*>(mcJet->TrackAt(k, mcArray));
984  if(!mcConstituent)
985  continue;
986  if(mcConstituent->GetLabel() == constituent->GetLabel())
987  {
988  foundInMC = kTRUE;
989  break;
990  }
991  }
992 
993  if(foundInMC)
994  {
995  tracksFromMC++;
996  ptFromMC += constituent->Pt();
997  }
998  tracksTotal++;
999  ptTotal += constituent->Pt();
1000 
1001  }
1002  }
1003  else // Calculate MC track percentage from all tracks in MC
1004  {
1005  for(Int_t j = 0; j < jet->GetNumberOfTracks(); j++)
1006  {
1007  AliVParticle* constituent = static_cast<AliVParticle*>(jet->TrackAt(j, fTracksCont->GetArray()));
1008  if(!constituent)
1009  continue;
1010 
1011  // Plots on MC percentage in jet
1012  if(constituent->GetLabel() > 10000)
1013  {
1014  tracksFromMC++;
1015  ptFromMC += constituent->Pt();
1016  }
1017  tracksTotal++;
1018  ptTotal += constituent->Pt();
1019  }
1020  }
1021 
1022 
1023  trackRatio = 0.;
1024  if(tracksTotal)
1025  trackRatio = tracksFromMC/((Double_t)tracksTotal);
1026 
1028  ptTotal = jet->Pt() - fJetsCont->GetRhoVal()*jet->Area();
1029 
1030  ptRatio = 0.;
1031  if(ptTotal)
1032  ptRatio = ptFromMC/ptTotal;
1033 
1034 }
1035 
1036 //________________________________________________________________________
1038 {
1039  Double_t leadingVetoJetPt = -999.;
1040  AliEmcalJet* leadingVetoJet = 0;
1041  // Search for the 'leading' overlapping jet in the veto sample
1042  if(fJetVetoArray && jet)
1043  {
1044  for(Int_t j=0; j<fJetVetoArray->GetEntries(); j++)
1045  {
1046  UInt_t dummy = 0;
1047  AliEmcalJet* vetoJet = static_cast<AliEmcalJet*>(fJetVetoArray->At(j));
1048  // Check if veto jet is accepted
1049  if(!fJetsCont->AcceptJet(vetoJet , dummy))
1050  continue;
1051 
1052  // Check matching distance
1053  Double_t vetoPt = vetoJet->Pt() - vetoJet->Area()*fJetsCont->GetRhoVal();
1054  Double_t deltaEta = (jet->Eta()-vetoJet->Eta());
1055  Double_t deltaPhi = TMath::Min(TMath::Abs(jet->Phi()-vetoJet->Phi()),TMath::TwoPi() - TMath::Abs(jet->Phi()-vetoJet->Phi()));
1056  Double_t deltaR = TMath::Sqrt(deltaEta*deltaEta + deltaPhi*deltaPhi);
1057 
1058  if ((vetoPt > leadingVetoJetPt) && (deltaR <= fJetEmbeddingMaxDistance))
1059  {
1060  leadingVetoJetPt = vetoPt;
1061  leadingVetoJet = vetoJet;
1062  }
1063  }
1064  }
1065 
1066  return leadingVetoJet;
1067 }
1068 
1069 //________________________________________________________________________
1071 {
1072  Double_t leadingVetoJetPt = -999.;
1073  AliEmcalJet* leadingVetoJet = 0;
1074  // Search for the 'leading' jet in the veto sample
1075  if(fJetVetoArray)
1076  {
1077  for(Int_t j=0; j<fJetVetoArray->GetEntries(); j++)
1078  {
1079  UInt_t dummy = 0;
1080  AliEmcalJet* vetoJet = static_cast<AliEmcalJet*>(fJetVetoArray->At(j));
1081  // Check if veto jet is accepted
1082  if(!fJetsCont->AcceptJet(vetoJet , dummy))
1083  continue;
1084 
1085  // Check matching distance
1086  Double_t vetoPt = vetoJet->Pt() - vetoJet->Area()*fJetsCont->GetRhoVal();
1087  if (vetoPt > leadingVetoJetPt)
1088  {
1089  leadingVetoJetPt = vetoPt;
1090  leadingVetoJet = vetoJet;
1091  }
1092  }
1093  }
1094  return leadingVetoJet;
1095 }
1096 
1097 //________________________________________________________________________
1099 {
1100  // This is to use a full cone in phi even at the edges of phi (2pi -> 0) (0 -> 2pi)
1101  Double_t trackPhi = 0.0;
1102  if (track->Phi() > (TMath::TwoPi() - (radius-phi)))
1103  trackPhi = track->Phi() - TMath::TwoPi();
1104  else if (track->Phi() < (phi+radius - TMath::TwoPi()))
1105  trackPhi = track->Phi() + TMath::TwoPi();
1106  else
1107  trackPhi = track->Phi();
1108 
1109  if ( TMath::Abs(trackPhi-phi)*TMath::Abs(trackPhi-phi) + TMath::Abs(track->Eta()-eta)*TMath::Abs(track->Eta()-eta) <= radius*radius)
1110  return kTRUE;
1111 
1112  return kFALSE;
1113 }
1114 
1115 
1116 //________________________________________________________________________
1118 {
1119  // Calculate leading + subleading jet
1121  if(fJetOutputMode==6)
1123  else if(fJetOutputMode==4 || fJetOutputMode==5)
1124  GetMatchingJets();
1125 }
1126 
1127 //________________________________________________________________________
1128 void AliAnalysisTaskChargedJetsHadronCF::GetLeadingJets(const char* opt, AliEmcalJet*& jetLeading, AliEmcalJet*& jetSubLeading)
1129 {
1130  // Customized from AliJetContainer::GetLeadingJet()
1131  // Get the leading+subleading jet; if opt contains "rho" the sorting is according to pt-A*rho
1132 
1133  TString option(opt);
1134  option.ToLower();
1135 
1136  jetLeading = 0;
1137  jetSubLeading = 0;
1138 
1139  fJetsCont->ResetCurrentID();
1140  Double_t tmpLeadingPt = 0;
1141  Double_t tmpSubleadingPt = 0;
1142 
1143  if (option.Contains("rho")) {
1144  while (AliEmcalJet* jet = fJetsCont->GetNextAcceptJet()) {
1145  if ( (jet->Pt()-jet->Area()*fJetsCont->GetRhoVal()) > tmpLeadingPt )
1146  {
1147  jetSubLeading = jetLeading;
1148  jetLeading = jet;
1149  tmpSubleadingPt = tmpLeadingPt;
1150  tmpLeadingPt = jet->Pt()-jet->Area()*fJetsCont->GetRhoVal();
1151  }
1152  else if ( (jet->Pt()-jet->Area()*fJetsCont->GetRhoVal()) > tmpSubleadingPt )
1153  {
1154  jetSubLeading = jet;
1155  tmpSubleadingPt = jet->Pt()-jet->Area()*fJetsCont->GetRhoVal();
1156  }
1157  }
1158  }
1159  else {
1160  while (AliEmcalJet* jet = fJetsCont->GetNextAcceptJet()) {
1161  if ( (jet->Pt()) > tmpLeadingPt )
1162  {
1163  jetSubLeading = jetLeading;
1164  jetLeading = jet;
1165  tmpSubleadingPt = tmpLeadingPt;
1166  tmpLeadingPt = jet->Pt();
1167  }
1168  else if ( (jet->Pt()) > tmpSubleadingPt )
1169  {
1170  jetSubLeading = jet;
1171  tmpSubleadingPt = jet->Pt();
1172  }
1173  }
1174  }
1175 }
1176 
1177 //________________________________________________________________________
1178 void AliAnalysisTaskChargedJetsHadronCF::GetLeadingJetsInArray(TClonesArray* arr, const char* opt, AliEmcalJet*& jetLeading, AliEmcalJet*& jetSubLeading)
1179 {
1180  // Customized from AliJetContainer::GetLeadingJet()
1181  // Get the leading+subleading jet; if opt contains "rho" the sorting is according to pt-A*rho
1182 
1183  TString option(opt);
1184  option.ToLower();
1185 
1186  jetLeading = 0;
1187  jetSubLeading = 0;
1188  Double_t tmpLeadingPt = -999.;
1189  Double_t tmpSubleadingPt = -999.;
1190 
1191  for(Int_t i=0; i<arr->GetEntries(); i++)
1192  {
1193  AliEmcalJet* jet = static_cast<AliEmcalJet*>(arr->At(i));
1194  UInt_t dummy = 0;
1195  if(!fJetsCont->AcceptJet(jet , dummy))
1196  continue;
1197 
1198  Double_t jetPt = jet->Pt();
1199  if (option.Contains("rho"))
1200  jetPt -= jet->Area()*fJetsCont->GetRhoVal();
1201 
1202  if ( jetPt > tmpLeadingPt )
1203  {
1204  jetSubLeading = jetLeading;
1205  jetLeading = jet;
1206  tmpSubleadingPt = tmpLeadingPt;
1207  tmpLeadingPt = jetPt;
1208  }
1209  else if ( jetPt > tmpSubleadingPt )
1210  {
1211  jetSubLeading = jet;
1212  tmpSubleadingPt = jetPt;
1213  }
1214 
1215  }
1216 }
1217 
1218 
1219 //________________________________________________________________________
1221 {
1222  // Method for the correct logarithmic binning of histograms
1223  TAxis *axis = h->GetAxis(axisNumber);
1224  int bins = axis->GetNbins();
1225 
1226  Double_t from = axis->GetXmin();
1227  Double_t to = axis->GetXmax();
1228  Double_t *newBins = new Double_t[bins + 1];
1229 
1230  newBins[0] = from;
1231  Double_t factor = pow(to/from, 1./bins);
1232 
1233  for (int i = 1; i <= bins; i++) {
1234  newBins[i] = factor * newBins[i-1];
1235  }
1236  axis->Set(bins, newBins);
1237  delete [] newBins;
1238 }
1239 
1240 //________________________________________________________________________
1242 {
1243  std::vector<AliEmcalJet*>::iterator matchedJetFindResult = std::find(fMatchedJets.begin(), fMatchedJets.end(), jet);
1244  if(matchedJetFindResult == fMatchedJets.end())
1245  {
1246  AliError("Checked for a reference jet but it was not found. Check code.");
1247  return 0;
1248  }
1249 
1250  Int_t matchedJetIndex = (matchedJetFindResult - fMatchedJets.begin());
1251  AliEmcalJet* refJet = fMatchedJetsReference[matchedJetIndex];
1252  return refJet;
1253 }
1254 
1255 //________________________________________________________________________
1257 {
1258  TH1* tmpHist = static_cast<TH1*>(fOutput->FindObject(key));
1259  if(!tmpHist)
1260  {
1261  AliError(Form("Cannot find histogram <%s> ",key)) ;
1262  return;
1263  }
1264 
1265  tmpHist->Fill(x);
1266 }
1267 
1268 //________________________________________________________________________
1270 {
1271  TH1* tmpHist = static_cast<TH1*>(fOutput->FindObject(key));
1272  if(!tmpHist)
1273  {
1274  AliError(Form("Cannot find histogram <%s> ",key));
1275  return;
1276  }
1277 
1278  if (tmpHist->IsA()->GetBaseClass("TH1"))
1279  static_cast<TH1*>(tmpHist)->Fill(x,y); // Fill x with y
1280  else if (tmpHist->IsA()->GetBaseClass("TH2"))
1281  static_cast<TH2*>(tmpHist)->Fill(x,y); // Fill x,y with 1
1282 }
1283 
1284 //________________________________________________________________________
1286 {
1287  TH2* tmpHist = static_cast<TH2*>(fOutput->FindObject(key));
1288  if(!tmpHist)
1289  {
1290  AliError(Form("Cannot find histogram <%s> ",key));
1291  return;
1292  }
1293 
1294  tmpHist->Fill(x,y,add);
1295 }
1296 
1297 //________________________________________________________________________
1299 {
1300  TH3* tmpHist = static_cast<TH3*>(fOutput->FindObject(key));
1301  if(!tmpHist)
1302  {
1303  AliError(Form("Cannot find histogram <%s> ",key));
1304  return;
1305  }
1306 
1307  if(add)
1308  tmpHist->Fill(x,y,z,add);
1309  else
1310  tmpHist->Fill(x,y,z);
1311 }
1312 
1313 
1314 //________________________________________________________________________
1315 template <class T> T* AliAnalysisTaskChargedJetsHadronCF::AddHistogram1D(const char* name, const char* title, const char* options, Int_t xBins, Double_t xMin, Double_t xMax, const char* xTitle, const char* yTitle)
1316 {
1317  T* tmpHist = new T(name, title, xBins, xMin, xMax);
1318 
1319  tmpHist->GetXaxis()->SetTitle(xTitle);
1320  tmpHist->GetYaxis()->SetTitle(yTitle);
1321  tmpHist->SetOption(options);
1322  tmpHist->SetMarkerStyle(kFullCircle);
1323  tmpHist->Sumw2();
1324 
1325  fOutput->Add(tmpHist);
1326 
1327  return tmpHist;
1328 }
1329 
1330 //________________________________________________________________________
1331 template <class T> T* AliAnalysisTaskChargedJetsHadronCF::AddHistogram2D(const char* name, const char* title, const char* options, Int_t xBins, Double_t xMin, Double_t xMax, Int_t yBins, Double_t yMin, Double_t yMax, const char* xTitle, const char* yTitle, const char* zTitle)
1332 {
1333  T* tmpHist = new T(name, title, xBins, xMin, xMax, yBins, yMin, yMax);
1334  tmpHist->GetXaxis()->SetTitle(xTitle);
1335  tmpHist->GetYaxis()->SetTitle(yTitle);
1336  tmpHist->GetZaxis()->SetTitle(zTitle);
1337  tmpHist->SetOption(options);
1338  tmpHist->SetMarkerStyle(kFullCircle);
1339  tmpHist->Sumw2();
1340 
1341  fOutput->Add(tmpHist);
1342 
1343  return tmpHist;
1344 }
1345 
1346 //________________________________________________________________________
1347 template <class T> T* AliAnalysisTaskChargedJetsHadronCF::AddHistogram3D(const char* name, const char* title, const char* options, Int_t xBins, Double_t xMin, Double_t xMax, Int_t yBins, Double_t yMin, Double_t yMax, Int_t zBins, Double_t zMin, Double_t zMax, const char* xTitle, const char* yTitle, const char* zTitle)
1348 {
1349  T* tmpHist = new T(name, title, xBins, xMin, xMax, yBins, yMin, yMax, zBins, zMin, zMax);
1350  tmpHist->GetXaxis()->SetTitle(xTitle);
1351  tmpHist->GetYaxis()->SetTitle(yTitle);
1352  tmpHist->GetZaxis()->SetTitle(zTitle);
1353  tmpHist->SetOption(options);
1354  tmpHist->SetMarkerStyle(kFullCircle);
1355  tmpHist->Sumw2();
1356 
1357  fOutput->Add(tmpHist);
1358 
1359  return tmpHist;
1360 }
1361 
1362 //________________________________________________________________________
1364 {
1365  // Called once at the end of the analysis.
1366 }
1367 
Int_t fArrayIndex
array index that holds the output array index
Double_t fEventExtractionMaxJetPt
maximum jet pt of recorded events
Short_t Charge() const
Definition: AliEmcalJet.h:123
Double_t Area() const
Definition: AliEmcalJet.h:130
void ExecOnce()
Perform steps needed to initialize the analysis.
Double_t GetRhoVal() const
double Double_t
Definition: External.C:58
Bool_t fJetVetoJetByJet
If true, the jet veto will be applied on a jet-by-jet basis.
Double_t GetJetEtaMin() const
T * AddHistogram1D(const char *name="CustomHistogram", const char *title="NO_TITLE", const char *options="", Int_t xBins=100, Double_t xMin=0.0, Double_t xMax=20.0, const char *xTitle="x axis", const char *yTitle="y axis")
const char * title
Definition: MakeQAPdf.C:27
void GetLeadingJets(const char *opt, AliEmcalJet *&jetLeading, AliEmcalJet *&jetSubLeading)
Float_t GetPartonEta6() const
AliEmcalPythiaInfo * fPythiaInfo
!event parton info
AliJetContainer * GetJetContainer(Int_t i=0) const
Definition: External.C:244
Float_t GetPartonEta7() const
std::vector< AliChargedJetsHadronCFCuts > fJetEmbeddingCuts
Cuts used in jet embedding.
Double_t Eta() const
Definition: AliEmcalJet.h:121
Double_t Phi() const
Definition: AliEmcalJet.h:117
Container with name, TClonesArray and cuts for particles.
Support task for (charged) jet-hadron correlations.
AliBasicParticle * fTreeBufferTrack
! Tree of extracted jets (buffer)
Double_t GetJetEtaMax() const
Float_t GetPartonPhi7() const
virtual Bool_t AcceptJet(Int_t i, UInt_t &rejectionReason) const
Int_t fJetOutputMode
mode which jets are written to array (0: all accepted, 1: leading, 2: subleading, 3: leading+subleadi...
Double_t fEventExtractionMinJetPt
minimum jet pt of recorded events
Bool_t fUseMCConstituents
If true, tracks with labels >= 10000 will be processed.
std::vector< AliEmcalJet * > fMatchedJetsReference
Jets matched in an event (reference)
Double_t fEventExtractionPercentage
percentage of events that is recorded
std::vector< TClonesArray * > fJetsOutput
! vector of arrays of basic correlation particles attached to the event (jets)
AliEmcalJet * fSubleadingJet
! subleading jet (calculated event-by-event)
void FillHistogram3D(const char *key, Double_t x, Double_t y, Double_t z, Double_t add=0)
Float_t GetPartonPhi6() const
Int_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:160
Container class of cuts for AliAnalysisTaskChargedJetsHadronCF.
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:139
UShort_t T(UShort_t m, UShort_t t)
Definition: RingBits.C:60
void GetLeadingJetsInArray(TClonesArray *arr, const char *opt, AliEmcalJet *&jetLeading, AliEmcalJet *&jetSubLeading)
TClonesArray * fJetEmbeddingArray
! Array of generated jets imported into task (for embedding)
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
void GetTrackMCRatios(AliEmcalJet *jet, AliEmcalJet *mcJet, Double_t &trackRatio, Double_t &ptRatio)
AliParticleContainer * GetParticleContainer() const
Double_t fTrackExtractionPercentagePower
Extraction percentage for tracks.
Bool_t fJetEmbeddingCreatePtPlotPerCut
create TH3 per cut or only once
int Int_t
Definition: External.C:63
AliAODPid * fTreeBufferPID
! Tree of extracted jets (buffer)
unsigned int UInt_t
Definition: External.C:33
TString fJetParticleArrayName
Name of fJetsOutput array (if one uses only one)
Double_t GetLeadingHadronPt(const AliEmcalJet *jet) const
T * AddHistogram3D(const char *name="CustomHistogram", const char *title="NO_TITLE", const char *options="", Int_t xBins=100, Double_t xMin=0.0, Double_t xMax=20.0, Int_t yBins=100, Double_t yMin=0.0, Double_t yMax=20.0, Int_t zBins=100, Double_t zMin=0.0, Double_t zMax=20.0, const char *xTitle="x axis", const char *yTitle="y axis", const char *zTitle="z axis")
TString fJetEmbeddingArrayName
Name of array used to match jets.
Int_t fNumberOfCentralityBins
Number of centrality bins.
Int_t fTreeBufferPDG
! Tree of extracted jets (buffer)
virtual AliVParticle * GetLeadingParticle(const char *opt="")
std::vector< AliEmcalJet * > fMatchedJets
Jets matched in an event (embedded)
TString fJetEmbeddingTrackArrayName
Name of array used to match tracks of jets.
Bool_t fUseDataConstituents
If true, tracks with labels < 10000 will be processed.
Double_t fCent
!event centrality
TString fTrackParticleArrayName
Name of fTracksOutput array.
TClonesArray * fTracksOutput
! Array of basic correlation particles attached to the event (tracks)
AliEmcalJet * fLeadingJet
! leading jet (calculated event-by-event)
AliEmcalJet * GetNextAcceptJet()
void AddJetToOutputArray(AliEmcalJet *jet, Int_t arrayIndex, Int_t &jetsAlreadyInArray)
Double_t Pt() const
Definition: AliEmcalJet.h:109
Bool_t IsTrackInCone(AliVParticle *track, Double_t eta, Double_t phi, Double_t radius)
Bool_t fJetEmbeddingUsePerTrackMCPercentage
When cutting on MC percentage, calculate it per track and not for all MC tracks.
Int_t fAcceptedJets
temporary var that holds how many jets passed
TString fJetVetoArrayName
Name of array used for veto jets.
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.
Float_t GetJetRadius() const
AliEmcalList * fOutput
!output list
Int_t fJetEmbeddingNumMatchedJets
Number of matched leading jets that will be used.
AliEmcalJet * fInitialPartonMatchedJet1
! On PYTHIA data and fJetOutputMode=6, this holds the PDG code of the initial collisions that was mat...
AliEmcalJet * fInitialPartonMatchedJet2
! On PYTHIA data and fJetOutputMode=6, this holds the PDG code of the initial collisions that was mat...
Definition: External.C:220
TH1 * fHistEventRejection
!book keep reasons for rejecting event
void SetMakeGeneralHistograms(Bool_t g)
Base task in the EMCAL jet framework.
Bool_t IsCutFulfilled(Double_t pt, Double_t mcPt, Double_t cent, Double_t ptRatio, Double_t vetoPt)
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:51
Bool_t fJetEmbeddingUseBgrdForMCPercentage
When cutting on MC percentage, use bgrd. corr to calculate MC percentage.
Declaration of class AliEmcalPythiaInfo.
const char Option_t
Definition: External.C:48
void UserCreateOutputObjects()
Main initialization function on the worker.
virtual AliVParticle * GetNextAcceptParticle()
bool Bool_t
Definition: External.C:53
Double_t fJetEmbeddingMaxDistance
Max distance allowed to accept an embedded jet.
void FillHistogramsJets(AliEmcalJet *jet, const char *cutName)
Int_t fAcceptedJets
! number accepted jets (calculated event-by-event)
T * AddHistogram2D(const char *name="CustomHistogram", const char *title="NO_TITLE", const char *options="", Int_t xBins=100, Double_t xMin=0.0, Double_t xMax=20.0, Int_t yBins=100, Double_t yMin=0.0, Double_t yMax=20.0, const char *xTitle="x axis", const char *yTitle="y axis", const char *zTitle="z axis")
Definition: External.C:196
TClonesArray * fJetVetoArray
! Array of jets imported into task used for veto a matching/embedding
Int_t fAcceptedTracks
! number accepted tracks (calculated event-by-event)
Int_t fNumRandomConesPerEvent
Number of random cones thrown in one event.
Int_t fConstPtFilterBit
For const pt plot, filter bit.