AliPhysics  cdeda5a (cdeda5a)
 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 "AliPicoTrack.h"
42 #include "AliVParticle.h"
43 #include "TRandom3.h"
45 
47 
48 //________________________________________________________________________
50 {
51 // dummy destructor
52 }
53 
54 //________________________________________________________________________
56 {
57 // dummy destructor
58 }
59 
60 //________________________________________________________________________
62 {
63 // dummy destructor
64 }
65 
69 //________________________________________________________________________
72  fJetsCont(0),
73  fTracksCont(0),
74  fJetsTree(0),
75  fJetsTreeBuffer(0),
76  fExtractionPercentage(0),
77  fExtractionMinPt(0),
78  fExtractionMaxPt(0),
79  fEventExtractionPercentage(0),
80  fEventExtractionMinJetPt(0),
81  fEventExtractionMaxJetPt(0),
82  fHadronMatchingRadius(0.5),
83  fConstPtFilterBit(1024),
84  fNumberOfCentralityBins(10),
85  fJetsOutput(),
86  fTracksOutput(),
87  fJetParticleArrayName("JetsDPhiBasicParticles"),
88  fTrackParticleArrayName(""),
89  fJetEmbeddingArray(),
90  fJetEmbeddingArrayName(""),
91  fJetEmbeddingTrackArrayName(""),
92  fJetEmbeddingMaxDistance(0.3),
93  fJetEmbeddingNumMatchedJets(2),
94  fJetEmbeddingUsePerTrackMCPercentage(kTRUE),
95  fJetEmbeddingUseBgrdForMCPercentage(kFALSE),
96  fJetEmbeddingCreatePtPlotPerCut(kFALSE),
97  fJetEmbeddingCuts(),
98  fJetVetoArray(),
99  fJetVetoArrayName(""),
100  fJetVetoJetByJet(1),
101  fMatchedJets(),
102  fRandom(0),
103  fJetOutputMode(0),
104  fPythiaExtractionMode(0),
105  fPythiaExtractionUseHadronMatching(kFALSE),
106  fLeadingJet(0),
107  fSubleadingJet(0),
108  fInitialPartonMatchedJet1(0),
109  fInitialPartonMatchedJet2(0),
110  fAcceptedJets(0),
111  fAcceptedTracks(0)
112 {
113  // Default constructor.
114  SetMakeGeneralHistograms(kTRUE);
115  fRandom = new TRandom3(0);
116 }
117 
118 
119 //________________________________________________________________________
121  AliAnalysisTaskEmcalJet(name, kTRUE),
122  fJetsCont(0),
123  fTracksCont(0),
124  fJetsTree(0),
125  fJetsTreeBuffer(0),
126  fExtractionPercentage(0),
127  fExtractionMinPt(0),
128  fExtractionMaxPt(0),
129  fEventExtractionPercentage(0),
130  fEventExtractionMinJetPt(0),
131  fEventExtractionMaxJetPt(0),
132  fHadronMatchingRadius(0.5),
133  fConstPtFilterBit(1024),
134  fNumberOfCentralityBins(10),
135  fJetsOutput(),
136  fTracksOutput(),
137  fJetParticleArrayName("JetsDPhiBasicParticles"),
138  fTrackParticleArrayName(""),
139  fJetEmbeddingArray(),
140  fJetEmbeddingArrayName(""),
141  fJetEmbeddingTrackArrayName(""),
142  fJetEmbeddingMaxDistance(0.3),
143  fJetEmbeddingNumMatchedJets(2),
144  fJetEmbeddingUsePerTrackMCPercentage(kTRUE),
145  fJetEmbeddingUseBgrdForMCPercentage(kFALSE),
146  fJetEmbeddingCreatePtPlotPerCut(kFALSE),
147  fJetEmbeddingCuts(),
148  fJetVetoArray(),
149  fJetVetoArrayName(""),
150  fJetVetoJetByJet(1),
151  fMatchedJets(),
152  fRandom(0),
153  fJetOutputMode(0),
154  fPythiaExtractionMode(0),
155  fPythiaExtractionUseHadronMatching(kFALSE),
156  fLeadingJet(0),
157  fSubleadingJet(0),
158  fInitialPartonMatchedJet1(0),
159  fInitialPartonMatchedJet2(0),
160  fAcceptedJets(0),
161  fAcceptedTracks(0)
162 {
163  // Constructor
165  fRandom = new TRandom3(0);
166 }
167 
168 //________________________________________________________________________
170 {
171  // Destructor.
172 }
173 
174 //________________________________________________________________________
176 {
178 
179  // ### Basic container settings
181  if(fJetsCont) { //get particles connected to jets
182  fJetsCont->PrintCuts();
184  } else { //no jets, just analysis tracks
186  }
187  if(fTracksCont) fTracksCont->SetClassName("AliVTrack");
188 
189  // ### Create all histograms
190  fHistEventRejection->GetXaxis()->SetBinLabel(15,"JetVeto");
191 
192  // Track QA plots
193  AddHistogram2D<TH2D>("hTrackPt", "Tracks p_{T} distribution", "", 300, 0., 300., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)", "Centrality", "dN^{Tracks}/dp_{T}");
194  AddHistogram2D<TH2D>("hTrackPhi", "Track angular distribution in #phi", "LEGO2", 180, 0., 2*TMath::Pi(), fNumberOfCentralityBins, 0, 100, "#phi", "Centrality", "dN^{Tracks}/(d#phi)");
195  AddHistogram2D<TH2D>("hTrackEta", "Track angular distribution in #eta", "LEGO2", 100, -2.5, 2.5, fNumberOfCentralityBins, 0, 100, "#eta", "Centrality", "dN^{Tracks}/(d#eta)");
196  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");
197 
198  AddHistogram2D<TH2D>("hLeadingTrackPt", "Leading tracks p_{T} distribution", "", 300, 0., 300., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)", "Centrality", "dN^{Tracks}/dp_{T}");
199  AddHistogram2D<TH2D>("hLeadingTrackPhi", "Leading tracks angular distribution in #phi", "LEGO2", 180, 0., 2*TMath::Pi(), fNumberOfCentralityBins, 0, 100, "#phi", "Centrality", "dN^{Tracks}/(d#phi)");
200  AddHistogram2D<TH2D>("hLeadingTrackEta", "Leading tracks angular distribution in #eta", "LEGO2", 100, -2.5, 2.5, fNumberOfCentralityBins, 0, 100, "#eta", "Centrality", "dN^{Tracks}/(d#eta)");
201  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");
202 
203  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})");
204  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})");
205 
206  // Create plots for each embedding cut
207  for(Int_t i = -1; i<static_cast<Int_t>(fJetEmbeddingCuts.size()); i++)
208  {
209  const char* appendix = "";
210  if(i>-1)
211  {
213  appendix = Form("_%s", currentCut.fCutName.Data());
214 
215  // Don't double-add cuts
216  if( static_cast<TH1*>(fOutput->FindObject(Form("hJetPtRaw%s", appendix))) )
217  continue;
218  }
219  // Jet QA plots
220  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}");
221  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}");
222  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");
223  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");
224  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}");
225  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}");
226  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");
227  AddHistogram2D<TH2D>(Form("hJetArea%s", appendix), "Jet area", "LEGO2", 200, 0., 2., fNumberOfCentralityBins, 0, 100, "Jet A", "Centrality", "dN^{Jets}/dA");
228  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}");
229  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}");
230 
231  // Leading/subleading ...
232  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}");
233  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}");
234  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");
235  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");
236  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}");
237  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}");
238  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");
239  AddHistogram2D<TH2D>(Form("hLeadingJetArea%s", appendix), "Jet area", "LEGO2", 200, 0., 2., fNumberOfCentralityBins, 0, 100, "Jet A", "Centrality", "dN^{Jets}/dA");
240  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}");
241  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}");
242 
243  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}");
244  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}");
245  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");
246  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");
247  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}");
248  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}");
249  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");
250  AddHistogram2D<TH2D>(Form("hSubleadingJetArea%s", appendix), "Jet area", "LEGO2", 200, 0., 2., fNumberOfCentralityBins, 0, 100, "Jet A", "Centrality", "dN^{Jets}/dA");
251  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}");
252  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}");
253 
254  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}");
255  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}");
256  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}");
257  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}");
258  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}");
259 
260  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}");
261  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}");
262  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}");
263  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}");
264  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}");
265 
266  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}");
267  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}");
268  }
269 
270  // Embedding plots
271  if(fJetOutputMode == 4 || fJetOutputMode == 5)
272  {
273  Double_t maxRatio = 1.;
275  maxRatio = 2.;
276 
277  for(Int_t i = -1; i<static_cast<Int_t>(fJetEmbeddingCuts.size()); i++)
278  {
279  const char* appendix = "";
280  if(i>-1)
281  {
283  appendix = Form("_%s", currentCut.fCutName.Data());
284 
285  // Don't double-add cuts
286  if( static_cast<TH1*>(fOutput->FindObject(Form("hEmbeddingDeltaR%s", appendix))) )
287  continue;
288  }
289  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");
290  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");
291  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");
292  AddHistogram1D<TH1D>(Form("hEmbeddingJetPt%s", appendix), "Embedded jets p_{T} distribution", "", 200, -50., 150., "p_{T, jet} (GeV/c)", "dN/dp_{T}");
293  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");
294 
296  {
297  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");
298  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");
299  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");
300  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");
301  }
302  else
303  {
304  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");
305  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");
306  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");
307  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");
308  }
309  }
310 
312  {
313  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");
314  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");
315  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");
316  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");
317  }
318  }
319 
320  // Random cone plots, background, ...
321  AddHistogram2D<TH2D>("hRandomConePt", "Random cone p_{T} distribution", "", 400, -100., 300., fNumberOfCentralityBins, 0, 100, "p_{T, cone} (GeV/c)", "Centrality", "dN^{Tracks}/dp_{T}");
322  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}");
323  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}");
324  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}");
325 
326  AddHistogram2D<TH2D>("hTrackCount", "Number of tracks in acceptance vs. centrality", "LEGO2", 500, 0., 5000., fNumberOfCentralityBins, 0, 100, "N tracks","Centrality", "dN^{Events}/dN^{Tracks}");
327  AddHistogram2D<TH2D>("hJetCount", "Number of jets in acceptance vs. centrality", "LEGO2", 100, 0., 100., fNumberOfCentralityBins, 0, 100, "N Jets","Centrality", "dN^{Events}/dN^{Jets}");
328  AddHistogram2D<TH2D>("hBackgroundPt", "Background p_{T} distribution", "", 150, 0., 150., fNumberOfCentralityBins, 0, 100, "Background p_{T} (GeV/c)", "Centrality", "dN^{Events}/dp_{T}");
329 
330 
331  for(Int_t i = -1; i<static_cast<Int_t>(fJetEmbeddingCuts.size()); i++)
332  {
333  const char* appendix = "";
334  if(i>-1)
335  {
337  appendix = Form("_%s", currentCut.fCutName.Data());
338 
339  // Don't double-add cuts
340  if( static_cast<TH1*>(fOutput->FindObject(Form("hBackgroundPtJetPt_Cent0_100%s", appendix))) )
341  continue;
342  }
343  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}");
344  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}");
345  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}");
346  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}");
347  }
348 
349  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
350 }
351 
352 
353 //________________________________________________________________________
355 
357 
358 
359  // ### Add the jets as basic correlation particles to the event
360  // This output object carries all accepted jets
361  if (!(fInputEvent->FindListObject(Form("%s", fJetParticleArrayName.Data()))))
362  {
363  fJetsOutput.push_back(new TClonesArray("AliPicoTrack"));
364  fJetsOutput.at(0)->SetName(fJetParticleArrayName.Data());
365  fInputEvent->AddObject(fJetsOutput.at(0));
366  }
367  else
368  AliFatal(Form("%s: Object with name %s already in event!", GetName(), Form("%s", fJetParticleArrayName.Data())));
369 
370  // These output objects carry all jets that pass certain cuts
371  if( (fJetOutputMode==4 || fJetOutputMode==5) )
372  {
373  // before, check if all given output names are OK
374  for(Int_t i = 0; i<fJetEmbeddingCuts.size(); i++)
375  if (fInputEvent->FindListObject(Form("%s", fJetEmbeddingCuts.at(i).fOutputName.Data())))
376  AliFatal(Form("%s: Object with name %s already in event!", GetName(), Form("%s", fJetEmbeddingCuts.at(i).fOutputName.Data())));
377 
378  for(Int_t i = 0; i<fJetEmbeddingCuts.size(); i++)
379  {
380  // If the cut demands a new output stream, add it
381  if (!fInputEvent->FindListObject(Form("%s", fJetEmbeddingCuts.at(i).fOutputName.Data())))
382  {
383  fJetsOutput.push_back(new TClonesArray("AliPicoTrack"));
384  fJetsOutput.at(fJetsOutput.size()-1)->SetName(fJetEmbeddingCuts.at(i).fOutputName.Data());
385  fInputEvent->AddObject(fJetsOutput.at(fJetsOutput.size()-1));
386  fJetEmbeddingCuts.at(i).fArrayIndex = fJetsOutput.size()-1;
387 
388  // Set the array indices for all cuts that use this output stream
389  for(Int_t j = 0; j<fJetEmbeddingCuts.size(); j++)
390  {
391  if(fJetEmbeddingCuts.at(j).fArrayIndex != -1)
392  continue;
393  if(fJetEmbeddingCuts.at(j).fOutputName != fJetEmbeddingCuts.at(i).fOutputName)
394  continue;
395  fJetEmbeddingCuts.at(j).fArrayIndex = fJetEmbeddingCuts.at(i).fArrayIndex;
396  }
397  }
398  }
399  }
400 
401  // ##############################################
402  // ##############################################
403 
404  // ### Add the tracks as basic correlation particles to the event (optional)
405  if(fTrackParticleArrayName != "")
406  {
407  if (!(fInputEvent->FindListObject(Form("%s", fTrackParticleArrayName.Data()))))
408  {
409  fTracksOutput = new TClonesArray("AliPicoTrack");
410  fTracksOutput->SetName(fTrackParticleArrayName.Data());
411  fInputEvent->AddObject(fTracksOutput);
412  }
413  else
414  AliFatal(Form("%s: Object with name %s already in event!", GetName(), Form("%s", fTrackParticleArrayName.Data())));
415  }
416 
417  // ### Import jets for embedding (optional)
418  if(fJetEmbeddingArrayName != "")
419  {
420  fJetEmbeddingArray = static_cast<TClonesArray*>(InputEvent()->FindListObject(Form("%s", fJetEmbeddingArrayName.Data())));
421  if(!fJetEmbeddingArray)
422  AliFatal(Form("Importing jets for embedding failed! Array '%s' not found!", fJetEmbeddingArrayName.Data()));
423  }
424  else if(fJetOutputMode==4 || fJetOutputMode==5)
425  AliFatal(Form("fJetEmbeddingArrayName must be set in jet output mode 4 or 5."));
426 
427  // ### Import veto jets for matching (optional)
428  if(fJetVetoArrayName != "")
429  {
430  fJetVetoArray = static_cast<TClonesArray*>(InputEvent()->FindListObject(Form("%s", fJetVetoArrayName.Data())));
431  if(!fJetVetoArray)
432  AliFatal(Form("Importing jets for veto failed! Array '%s' not found!", fJetVetoArrayName.Data()));
433  }
434 
435  // ### Jets tree (optional)
437  {
438  fJetsTree = new TTree("ExtractedJets", "ExtractedJets");
439  fJetsTree->Branch("Jets", "AliBasicJet", &fJetsTreeBuffer, 1000);
440  fOutput->Add(fJetsTree);
441  }
442 
443 
444 }
445 
446 //________________________________________________________________________
448 {
449  if(fJetOutputMode==3) // output leading&subleading jet
450  {
451  if((jet!=fLeadingJet) && (jet!=fSubleadingJet))
452  return kFALSE;
453  }
454  else if(fJetOutputMode==1) // output the leading jet
455  {
456  if(jet!=fLeadingJet)
457  return kFALSE;
458  }
459  else if(fJetOutputMode==2) // output the subleading jet
460  {
461  if(jet!=fSubleadingJet)
462  return kFALSE;
463  }
464  else if(fJetOutputMode==6)
465  {
467  return kFALSE;
468  }
469 
470  if(fJetOutputMode==4) // matching jets only
471  return (std::find(fMatchedJets.begin(), fMatchedJets.end(), jet) != fMatchedJets.end());
472  else if(fJetOutputMode==5) // non-matching jets only
473  return (std::find(fMatchedJets.begin(), fMatchedJets.end(), jet) == fMatchedJets.end());
474 
475  return kTRUE;
476 }
477 
478 //________________________________________________________________________
480 {
481  TString appendix("");
482  if(cutName)
483  appendix = Form("_%s", cutName);
484 
485  // All jets
486  FillHistogram(Form("hJetPtRaw%s", appendix.Data()), jet->Pt(), fCent);
487  FillHistogram(Form("hJetPt%s", appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), fCent);
488  FillHistogram(Form("hJetPhi%s", appendix.Data()), jet->Phi(), fCent);
489  FillHistogram(Form("hJetEta%s", appendix.Data()), jet->Eta(), fCent);
490  FillHistogram(Form("hJetEtaPt%s", appendix.Data()), jet->Eta(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
491  FillHistogram(Form("hJetPhiPt%s", appendix.Data()), jet->Phi(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
492  FillHistogram(Form("hJetPhiEta%s", appendix.Data()), jet->Phi(), jet->Eta());
493  FillHistogram(Form("hJetArea%s", appendix.Data()), jet->Area(), fCent);
494  FillHistogram(Form("hJetAreaPt%s", appendix.Data()), jet->Area(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
495  FillHistogram(Form("hJetPtLeadingHadron%s", appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), fJetsCont->GetLeadingHadronPt(jet));
496 
497  FillHistogram(Form("hBackgroundPtJetPt_Cent0_100%s", appendix.Data()), fJetsCont->GetRhoVal(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
498  if( (fCent >= 0) && (fCent < 10) )
499  FillHistogram(Form("hBackgroundPtJetPt_Cent0_10%s", appendix.Data()), fJetsCont->GetRhoVal(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
500 
501  // Leading jet plots
502  if(jet==fLeadingJet)
503  {
504  FillHistogram(Form("hLeadingJetPtRaw%s", appendix.Data()), jet->Pt(), fCent);
505  FillHistogram(Form("hLeadingJetPt%s", appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), fCent);
506  FillHistogram(Form("hLeadingJetPhi%s", appendix.Data()), jet->Phi(), fCent);
507  FillHistogram(Form("hLeadingJetEta%s", appendix.Data()), jet->Eta(), fCent);
508  FillHistogram(Form("hLeadingJetEtaPt%s", appendix.Data()), jet->Eta(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
509  FillHistogram(Form("hLeadingJetPhiPt%s", appendix.Data()), jet->Phi(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
510  FillHistogram(Form("hLeadingJetPhiEta%s", appendix.Data()), jet->Phi(), jet->Eta());
511  FillHistogram(Form("hLeadingJetArea%s", appendix.Data()), jet->Area(), fCent);
512  FillHistogram(Form("hLeadingJetAreaPt%s", appendix.Data()), jet->Area(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
513  FillHistogram(Form("hLeadingJetPtLeadingHadron%s", appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), fJetsCont->GetLeadingHadronPt(jet));
514  }
515 
516  // Subleading jet plot
517  else if(jet==fSubleadingJet)
518  {
519  FillHistogram(Form("hSubleadingJetPtRaw%s", appendix.Data()), jet->Pt(), fCent);
520  FillHistogram(Form("hSubleadingJetPt%s", appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), fCent);
521  FillHistogram(Form("hSubleadingJetPhi%s", appendix.Data()), jet->Phi(), fCent);
522  FillHistogram(Form("hSubleadingJetEta%s", appendix.Data()), jet->Eta(), fCent);
523  FillHistogram(Form("hSubleadingJetEtaPt%s", appendix.Data()), jet->Eta(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
524  FillHistogram(Form("hSubleadingJetPhiPt%s", appendix.Data()), jet->Phi(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
525  FillHistogram(Form("hSubleadingJetPhiEta%s", appendix.Data()), jet->Phi(), jet->Eta());
526  FillHistogram(Form("hSubleadingJetArea%s", appendix.Data()), jet->Area(), fCent);
527  FillHistogram(Form("hSubleadingJetAreaPt%s", appendix.Data()), jet->Area(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
528  FillHistogram(Form("hSubleadingJetPtLeadingHadron%s", appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), fJetsCont->GetLeadingHadronPt(jet));
529  }
530 
531  // ####### Jet constituent plots
532  for(Int_t i = 0; i < jet->GetNumberOfTracks(); i++)
533  {
534  AliVParticle* constituent = static_cast<AliVParticle*>(jet->TrackAt(i, fTracksCont->GetArray()));
535  if(!constituent)
536  continue;
537 
538  Bool_t filterConditionFulfilled = kFALSE;
539  AliAODTrack* aodTrack = static_cast<AliAODTrack*>(constituent);
540  if (aodTrack)
541  filterConditionFulfilled = aodTrack->TestFilterBit(fConstPtFilterBit);
542 
543  // Fill jet constituent plots
544  FillHistogram(Form("hJetConstituentPt_Cent0_100%s", appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), constituent->Pt());
545  if(filterConditionFulfilled)
546  FillHistogram(Form("hJetConstituentPt_Cent0_100_FilterBit%i%s", fConstPtFilterBit, appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), constituent->Pt());
547  if( (fCent >= 0) && (fCent < 10) )
548  {
549  FillHistogram(Form("hJetConstituentPt_Cent0_10%s", appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), constituent->Pt());
550  if(filterConditionFulfilled)
551  FillHistogram(Form("hJetConstituentPt_Cent0_10_FilterBit%i%s", fConstPtFilterBit, appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), constituent->Pt());
552  }
553  else if( (fCent >= 10) && (fCent < 30) )
554  {
555  FillHistogram(Form("hJetConstituentPt_Cent10_30%s", appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), constituent->Pt());
556  if(filterConditionFulfilled)
557  FillHistogram(Form("hJetConstituentPt_Cent10_30_FilterBit%i%s", fConstPtFilterBit, appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), constituent->Pt());
558  }
559  else if( (fCent >= 30) && (fCent < 50) )
560  {
561  FillHistogram(Form("hJetConstituentPt_Cent30_50%s", appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), constituent->Pt());
562  if(filterConditionFulfilled)
563  FillHistogram(Form("hJetConstituentPt_Cent30_50_FilterBit%i%s", fConstPtFilterBit, appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), constituent->Pt());
564  }
565  else if( (fCent >= 50) && (fCent < 90) )
566  {
567  FillHistogram(Form("hJetConstituentPt_Cent50_90%s", appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), constituent->Pt());
568  if(filterConditionFulfilled)
569  FillHistogram(Form("hJetConstituentPt_Cent50_90_FilterBit%i%s", fConstPtFilterBit, appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), constituent->Pt());
570  }
571  }
572 
573  FillHistogram(Form("hJetConstituentCount_Cent0_100%s", appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), jet->GetNumberOfTracks());
574  if( (fCent >= 0) && (fCent < 10) )
575  FillHistogram(Form("hJetConstituentCount_Cent0_10%s", appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), jet->GetNumberOfTracks());
576 
577  FillHistogram(Form("hBackgroundPtConstCount_Cent0_100%s", appendix.Data()), fJetsCont->GetRhoVal(), jet->GetNumberOfTracks());
578  if( (fCent >= 0) && (fCent < 10) )
579  FillHistogram(Form("hBackgroundPtConstCount_Cent0_10%s", appendix.Data()), fJetsCont->GetRhoVal(), jet->GetNumberOfTracks());
580 
581  // ####### Embedding plots
582  if( (fJetOutputMode == 4) || (fJetOutputMode == 5))
583  {
584  AliEmcalJet* refJet = GetReferenceJet(jet);
585  Double_t deltaEta = (jet->Eta()-refJet->Eta());
586  Double_t deltaPhi = TMath::Min(TMath::Abs(jet->Phi()-refJet->Phi()),TMath::TwoPi() - TMath::Abs(jet->Phi()-refJet->Phi()));
587  if(jet->Phi() < refJet->Phi())
588  deltaPhi = -deltaPhi;
589 
590  Double_t deltaR = TMath::Sqrt(deltaEta*deltaEta + deltaPhi*deltaPhi);
591  FillHistogram(Form("hEmbeddingDeltaR%s", appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), deltaR);
592  FillHistogram(Form("hEmbeddingDeltaEta%s", appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), deltaPhi);
593  FillHistogram(Form("hEmbeddingDeltaPhi%s", appendix.Data()), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), deltaEta);
594  FillHistogram(Form("hEmbeddingJetPt%s", appendix.Data()), refJet->Pt());
595  FillHistogram(Form("hEmbeddingJetPhiEta%s", appendix.Data()), refJet->Phi(), refJet->Eta());
596 
597  // Only create 3D plots for each cut on demand
598  Double_t trackRatio = 0.;
599  Double_t ptRatio = 0.;
600  GetTrackMCRatios(jet, refJet, trackRatio, ptRatio);
601 
602  if(fCent >= 0 && fCent < 10)
603  {
604  if((appendix == "") || fJetEmbeddingCreatePtPlotPerCut)
605  FillHistogram3D(Form("hEmbeddingPtCorr010%s", appendix.Data()), refJet->Pt(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), ptRatio);
606  if(ptRatio >= 0.2)
607  FillHistogram(Form("hEmbeddingPtCorr010_Above20%s", appendix.Data()), refJet->Pt(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
608  }
609  else if (fCent >= 10 && fCent < 30)
610  {
611  if((appendix == "") || fJetEmbeddingCreatePtPlotPerCut)
612  FillHistogram3D(Form("hEmbeddingPtCorr1030%s", appendix.Data()), refJet->Pt(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), ptRatio);
613  if(ptRatio >= 0.2)
614  FillHistogram(Form("hEmbeddingPtCorr1030_Above20%s", appendix.Data()), refJet->Pt(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
615  }
616  else if (fCent >= 30 && fCent < 50)
617  {
618  if((appendix == "") || fJetEmbeddingCreatePtPlotPerCut)
619  FillHistogram3D(Form("hEmbeddingPtCorr3050%s", appendix.Data()), refJet->Pt(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), ptRatio);
620  if(ptRatio >= 0.2)
621  FillHistogram(Form("hEmbeddingPtCorr3050_Above20%s", appendix.Data()), refJet->Pt(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
622  }
623  else if (fCent >= 50 && fCent < 90)
624  {
625  if((appendix == "") || fJetEmbeddingCreatePtPlotPerCut)
626  FillHistogram3D(Form("hEmbeddingPtCorr5090%s", appendix.Data()), refJet->Pt(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), ptRatio);
627  if(ptRatio >= 0.2)
628  FillHistogram(Form("hEmbeddingPtCorr5090_Above20%s", appendix.Data()), refJet->Pt(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
629  }
630  }
631 }
632 
633 //________________________________________________________________________
635 {
636  FillHistogram("hTrackPt", track->Pt(), fCent);
637  FillHistogram("hTrackPhi", track->Phi(), fCent);
638  FillHistogram("hTrackEta", track->Eta(), fCent);
639  FillHistogram("hTrackEtaPt", track->Eta(), track->Pt());
640  FillHistogram("hTrackPhiPt", track->Phi(), track->Pt());
641  FillHistogram("hTrackPhiEta", track->Phi(), track->Eta());
642 }
643 
644 //________________________________________________________________________
646 {
647  Double_t tmpPt = jet->Pt() - fJetsCont->GetRhoVal()*jet->Area();
648  new ((*fJetsOutput.at(arrayIndex))[jetsAlreadyInArray]) AliPicoTrack(tmpPt, jet->Eta(), jet->Phi(), jet->Charge(), 0, 0);
649  jetsAlreadyInArray++;
650 }
651 
652 //________________________________________________________________________
654 {
655  if(fTrackParticleArrayName != "")
656  {
657  new ((*fTracksOutput)[fAcceptedTracks]) AliPicoTrack(track->Pt(), track->Eta(), track->Phi(), track->Charge(), 0, 0); // only Pt,Eta,Phi are interesting for correlations;
658  fAcceptedTracks++;
659  }
660 }
661 
662 //________________________________________________________________________
664 {
665  // Check pT threshold
666  if( ((jet->Pt()-jet->Area()*fJetsCont->GetRhoVal()) < fExtractionMinPt) || ((jet->Pt()-jet->Area()*fJetsCont->GetRhoVal()) >= fExtractionMaxPt) )
667  return;
668 
669  AliVHeader* eventIDHeader = InputEvent()->GetHeader();
670  Long64_t eventID = 0;
671  if(eventIDHeader)
672  eventID = eventIDHeader->GetEventIdAsLong();
673 
674  // if only the two initial collision partons will be added, get info on them
675  Int_t matchedIC = 0;
676  Int_t matchedHadron = 0;
677  if(fJetOutputMode==6)
678  {
679  // Get type of jet
680  CalculateJetType(jet, matchedIC, matchedHadron);
681  Int_t partid = matchedIC;
683  partid = matchedHadron;
684  // If fPythiaExtractionMode is set, only extract certain jets
685  if( (fPythiaExtractionMode==1) && not (partid>=1 && partid<=6)) // all quark-jet extraction
686  return;
687  else if( (fPythiaExtractionMode==2) && not (partid==21)) // gluon-jet extraction
688  return;
689  else if( (fPythiaExtractionMode==3) && not (partid==0)) // extract only those w/o hadron matching
690  return;
691  else if( (fPythiaExtractionMode<0) && (fPythiaExtractionMode!=-partid) ) // custom type jet extraction by given a negative number
692  return;
693  }
694 
695  // Load vertex if possible
696  Double_t vtxX = 0;
697  Double_t vtxY = 0;
698  Double_t vtxZ = 0;
699  const AliVVertex* myVertex = InputEvent()->GetPrimaryVertex();
700  if(!myVertex && MCEvent())
701  myVertex = MCEvent()->GetPrimaryVertex();
702  if(myVertex)
703  {
704  vtxX = myVertex->GetX();
705  vtxY = myVertex->GetY();
706  vtxZ = myVertex->GetZ();
707  }
708 
709  // Discard jets statistically
710  if(fRandom->Rndm() >= fExtractionPercentage)
711  return;
712 
713  AliBasicJet basicJet(jet->Eta(), jet->Phi(), jet->Pt(), jet->Charge(), fJetsCont->GetJetRadius(), jet->Area(), matchedIC, matchedHadron, fJetsCont->GetRhoVal(), InputEvent()->GetMagneticField(), vtxX, vtxY, vtxZ, eventID, fCent);
714 
715  // Add constituents
716  for(Int_t i = 0; i < jet->GetNumberOfTracks(); i++)
717  {
718  AliVParticle* particle = static_cast<AliVParticle*>(jet->TrackAt(i, fTracksCont->GetArray()));
719  if(!particle) continue;
720 
721  AliAODTrack* aodtrack = static_cast<AliAODTrack*>(jet->TrackAt(i, fTracksCont->GetArray()));
722  Int_t constid = 9; // 9 mean unknown
723  if(fJetOutputMode==6) // MC
724  {
725  // Use same convention as PID in AODs
726  if(TMath::Abs(particle->PdgCode()) == 2212) // proton
727  constid = 4;
728  else if (TMath::Abs(particle->PdgCode()) == 211) // pion
729  constid = 2;
730  else if (TMath::Abs(particle->PdgCode()) == 321) // kaon
731  constid = 3;
732  else if (TMath::Abs(particle->PdgCode()) == 11) // electron
733  constid = 0;
734  else if (TMath::Abs(particle->PdgCode()) == 13) // muon
735  constid = 1;
736  }
737  else if (aodtrack) // data
738  constid = aodtrack->GetMostProbablePID();
739 
740  basicJet.AddJetConstituent(particle->Eta(), particle->Phi(), particle->Pt(), particle->Charge(), constid, particle->Xv(), particle->Yv(), particle->Zv());
741  }
742  if(std::find(fMatchedJets.begin(), fMatchedJets.end(), jet) != fMatchedJets.end()) // set the true pT from the matched jets (only possible in modes 4 & 7)
743  basicJet.SetTruePt(fMatchedJetsReference[std::find(fMatchedJets.begin(), fMatchedJets.end(), jet)-fMatchedJets.begin()]->Pt());
744 
745  fJetsTreeBuffer = &basicJet;
746  fJetsTree->Fill();
747 }
748 
749 //________________________________________________________________________
751 {
752  // Check jet pT threshold
754  return;
755 
756  // Discard jets statistically
757  if(fRandom->Rndm() >= fEventExtractionPercentage)
758  return;
759 
760  static Int_t numSavedEvents = 0;
761  numSavedEvents++;
762 
763 
764  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");
765  fTracksCont->ResetCurrentID();
766  while(AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle()))
767  FillHistogram(Form("Event%i", numSavedEvents), track->Phi(), track->Eta(), track->Pt());
768 
769 }
770 
771 //________________________________________________________________________
773 {
775 
776  // Jet veto:
777  // 1. jet-by-jet mode: veto is active if jets overlaps with a suitable jet
778  // 2. other mode: veto is active if suitable jet is in sample
779  AliEmcalJet* vetoJet = 0;
780  if(!fJetVetoJetByJet)
781  vetoJet = GetLeadingVetoJet();
782 
783  // ####### Jet loop
784  fAcceptedJets = 0;
785  for(Int_t i = 0; i<fJetEmbeddingCuts.size(); i++)
786  fJetEmbeddingCuts.at(i).fAcceptedJets = 0;
787  fJetsCont->ResetCurrentID();
788  while(AliEmcalJet *jet = fJetsCont->GetNextAcceptJet())
789  {
790  if(!IsJetSelected(jet))
791  continue;
792 
793  // Plots + output jets regardless of embedding cut
794  FillHistogramsJets(jet, 0);
796 
797  // Plots + output jets for each embedding cut
798  if(fJetVetoJetByJet)
799  vetoJet = GetVetoJet(jet);
800  for(Int_t i = 0; i<fJetEmbeddingCuts.size(); i++)
801  {
803  AliEmcalJet* refJet = GetReferenceJet(jet);
804  Double_t trackRatio = 0.;
805  Double_t ptRatio = 0.;
806  GetTrackMCRatios(jet, refJet, trackRatio, ptRatio);
807 
808  Double_t vetoJetPt = 0.;
809  if(vetoJet)
810  vetoJetPt = vetoJet->Pt();
811 
812  if(!currentCut.IsCutFulfilled(jet->Pt(), refJet->Pt(), fCent, ptRatio, vetoJetPt))
813  continue;
814 
815  FillHistogramsJets(jet, currentCut.fCutName.Data());
816  AddJetToOutputArray(jet, currentCut.fArrayIndex, currentCut.fAcceptedJets);
817  }
818 
819  // Add jet to output array
821  AddJetToTree(jet);
822  }
823 
824  // ####### Particle loop
825  // Throw random cone
826  Double_t tmpRandConeEta = fJetsCont->GetJetEtaMin() + fRandom->Rndm()*TMath::Abs(fJetsCont->GetJetEtaMax()-fJetsCont->GetJetEtaMin());
827  Double_t tmpRandConePhi = fRandom->Rndm()*TMath::TwoPi();
828  Double_t tmpRandConePt = 0; // to be determined
829  Double_t tmpRandConePt3GeV = 0; // to be determined
830 
831  fAcceptedTracks = 0;
832  fTracksCont->ResetCurrentID();
833  Int_t trackcount = 0;
834  while(AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle()))
835  {
836  // Track plots
837  FillHistogramsTracks(track);
838 
839  if(IsTrackInCone(track, tmpRandConeEta, tmpRandConePhi, fJetsCont->GetJetRadius()))
840  {
841  tmpRandConePt += track->Pt();
842  if (track->Pt() > 3.0)
843  tmpRandConePt3GeV += track->Pt();
844  }
845 
846  // Add track to output array
847  trackcount++;
848  AddTrackToOutputArray(track);
849  }
850 
851  // Add event to output tree
853  AddEventToTree();
854 
855  // ####### Event properties
856  FillHistogram("hRandomConePt", tmpRandConePt - fJetsCont->GetRhoVal()*fJetsCont->GetJetRadius()*fJetsCont->GetJetRadius()*TMath::Pi(), fCent);
857  FillHistogram("hRandomConePtCut3GeV", tmpRandConePt3GeV - fJetsCont->GetRhoVal()*fJetsCont->GetJetRadius()*fJetsCont->GetJetRadius()*TMath::Pi(), fCent);
858  FillHistogram("hRandomConeRawPt", tmpRandConePt, fCent);
859  FillHistogram("hRandomConeRawPtCut3GeV", tmpRandConePt3GeV, fCent);
860 
861  FillHistogram("hBackgroundPt", fJetsCont->GetRhoVal(), fCent);
862  FillHistogram("hJetCount", fAcceptedJets, fCent);
863  FillHistogram("hTrackCount", trackcount, fCent);
864  // NOTE: It is possible to use fTracksCont->GetLeadingParticle() since we do not apply additional track cuts
865  AliVTrack* leadTrack = static_cast<AliVTrack*>(fTracksCont->GetLeadingParticle());
866  if(leadTrack)
867  {
868  FillHistogram("hLeadingTrackPt", leadTrack->Pt(), fCent);
869  FillHistogram("hLeadingTrackPhi", leadTrack->Phi(), fCent);
870  FillHistogram("hLeadingTrackEta", leadTrack->Eta(), fCent);
871  FillHistogram("hLeadingTrackPhiEta", leadTrack->Phi(), leadTrack->Eta());
872  }
873 
874  return kTRUE;
875 }
876 
877 //########################################################################
878 // HELPERS
879 //########################################################################
880 
881 //________________________________________________________________________
883 {
884  if(!fPythiaInfo)
885  AliError("fPythiaInfo object not available. Is it activated with SetGeneratePythiaInfoObject()?");
886 
887  Double_t bestMatchDeltaR1 = 999.;
888  Double_t bestMatchDeltaR2 = 999.;
889 
890  fJetsCont->ResetCurrentID();
891  while(AliEmcalJet *jet = fJetsCont->GetNextAcceptJet())
892  {
893  // Check via geometrical matching if jet is connected to the initial collision
894  Double_t deltaEta1 = TMath::Abs(jet->Eta()-fPythiaInfo->GetPartonEta6());
895  Double_t deltaEta2 = TMath::Abs(jet->Eta()-fPythiaInfo->GetPartonEta7());
896  Double_t deltaPhi1 = TMath::Min(TMath::Abs(jet->Phi()-fPythiaInfo->GetPartonPhi6()),TMath::TwoPi() - TMath::Abs(jet->Phi()-fPythiaInfo->GetPartonPhi6()));
897  Double_t deltaPhi2 = TMath::Min(TMath::Abs(jet->Phi()-fPythiaInfo->GetPartonPhi7()),TMath::TwoPi() - TMath::Abs(jet->Phi()-fPythiaInfo->GetPartonPhi7()));
898 
899  Double_t deltaR1 = TMath::Sqrt(deltaEta1*deltaEta1 + deltaPhi1*deltaPhi1);
900  Double_t deltaR2 = TMath::Sqrt(deltaEta2*deltaEta2 + deltaPhi2*deltaPhi2);
901 
902  if(deltaR1 < bestMatchDeltaR1)
903  {
904  bestMatchDeltaR1 = deltaR1;
906  }
907  if(deltaR2 < bestMatchDeltaR2)
908  {
909  bestMatchDeltaR2 = deltaR2;
911  }
912  }
913 
914  if(bestMatchDeltaR1 > 0.3)
916  if(bestMatchDeltaR2 > 0.3)
918 }
919 
920 //________________________________________________________________________
922 {
923  fMatchedJets.clear();
924  fMatchedJetsReference.clear();
925 
926  AliEmcalJet* jetLeading = 0;
927  AliEmcalJet* jetSubLeading = 0;
928  GetLeadingJetsInArray(fJetEmbeddingArray, "rho", jetLeading, jetSubLeading);
929 
930  // ############ Search for all matches with leading/subleading jets
931  for(Int_t i=0; i<fJetEmbeddingNumMatchedJets; i++)
932  {
933  AliEmcalJet* probeJet = 0;
934  if(i==0)
935  probeJet = jetLeading;
936  else if(i==1)
937  probeJet = jetSubLeading;
938 
939  if(!probeJet)
940  continue;
941 
942  AliEmcalJet* matchedJet = 0;
943  AliEmcalJet* matchedJetReference = 0;
944  Double_t bestMatchDeltaR = 999.;
945  fJetsCont->ResetCurrentID();
946  // Loop over all embedded jets to find the best match
947  while(AliEmcalJet* embeddedJet = fJetsCont->GetNextAcceptJet())
948  {
949  Double_t deltaEta = (embeddedJet->Eta()-probeJet->Eta());
950  Double_t deltaPhi = TMath::Min(TMath::Abs(embeddedJet->Phi()-probeJet->Phi()),TMath::TwoPi() - TMath::Abs(embeddedJet->Phi()-probeJet->Phi()));
951  Double_t deltaR = TMath::Sqrt(deltaEta*deltaEta + deltaPhi*deltaPhi);
952 
953  // Cut jets too far away
954  if (deltaR > fJetEmbeddingMaxDistance)
955  continue;
956 
957  // Search for the best match
958  if(deltaR < bestMatchDeltaR)
959  {
960  bestMatchDeltaR = deltaR;
961  matchedJet = embeddedJet;
962  matchedJetReference = probeJet;
963  }
964  }
965  // Put matched jet to a list
966  if(matchedJet && matchedJetReference)
967  {
968  fMatchedJets.push_back(matchedJet);
969  fMatchedJetsReference.push_back(matchedJetReference);
970  }
971  }
972 }
973 
974 //________________________________________________________________________
976 {
977  typeIC = 0;
978  if(!fPythiaInfo)
979  AliError("fPythiaInfo object not available. Is it activated with SetGeneratePythiaInfoObject()?");
980  else if(jet==fInitialPartonMatchedJet1)
981  typeIC = fPythiaInfo->GetPartonFlag6();
982  else if (jet==fInitialPartonMatchedJet2)
983  typeIC = fPythiaInfo->GetPartonFlag7();
984 
985 
986  typeHM = 0;
987  AliStack* stack = MCEvent()->Stack();
988  // Go through the whole particle stack
989  for(Int_t i=0; i<stack->GetNtrack(); i++)
990  {
991  TParticle *part = stack->Particle(i);
992  if(!part) continue;
993 
994  // Check if particle is in a radius around the jet
995  Double_t rsquared = (part->Eta() - jet->Eta())*(part->Eta() - jet->Eta()) + (part->Phi() - jet->Phi())*(part->Phi() - jet->Phi());
997  {
998  continue;
999  }
1000 
1001  // Check if the particle has beauty, charm or strangeness
1002  // If it has beauty, we are done (exclusive definition)
1003  Int_t absPDG = TMath::Abs(part->GetPdgCode());
1004  // Particle has beauty
1005  if ((absPDG > 500 && absPDG < 600) || (absPDG > 5000 && absPDG < 6000))
1006  {
1007  typeHM = 5; // beauty
1008  break;
1009  }
1010  // Particle has charm
1011  else if ((absPDG > 400 && absPDG < 500) || (absPDG > 4000 && absPDG < 5000))
1012  typeHM = 4; // charm
1013  // Particle has strangeness: Only search for strangeness, if charm was not already found
1014  else if (typeHM != 4 && (absPDG > 300 && absPDG < 400) || (absPDG > 3000 && absPDG < 4000))
1015  typeHM = 3; // strange
1016  }
1017 }
1018 
1019 
1020 //________________________________________________________________________
1022 {
1023  Int_t tracksFromMC = 0;
1024  Int_t tracksTotal = 0;
1025  Double_t ptFromMC = 0;
1026  Double_t ptTotal = 0;
1027 
1028  if(fJetEmbeddingUsePerTrackMCPercentage) // Calculate MC track percentage from tracks of the matched jet
1029  {
1030  for(Int_t j = 0; j < jet->GetNumberOfTracks(); j++)
1031  {
1032  AliVParticle* constituent = static_cast<AliVParticle*>(jet->TrackAt(j, fTracksCont->GetArray()));
1033  if(!constituent)
1034  continue;
1035 
1036  TClonesArray* mcArray = static_cast<TClonesArray*>(InputEvent()->FindListObject(Form("%s", fJetEmbeddingTrackArrayName.Data())));
1037  Bool_t foundInMC = kFALSE;
1038  for(Int_t k=0; k<mcJet->GetNumberOfTracks(); k++)
1039  {
1040  AliVParticle* mcConstituent = static_cast<AliVParticle*>(mcJet->TrackAt(k, mcArray));
1041  if(!mcConstituent)
1042  continue;
1043  if(mcConstituent->GetLabel() == constituent->GetLabel())
1044  {
1045  foundInMC = kTRUE;
1046  break;
1047  }
1048  }
1049 
1050  if(foundInMC)
1051  {
1052  tracksFromMC++;
1053  ptFromMC += constituent->Pt();
1054  }
1055  tracksTotal++;
1056  ptTotal += constituent->Pt();
1057 
1058  }
1059  }
1060  else // Calculate MC track percentage from all tracks in MC
1061  {
1062  for(Int_t j = 0; j < jet->GetNumberOfTracks(); j++)
1063  {
1064  AliVParticle* constituent = static_cast<AliVParticle*>(jet->TrackAt(j, fTracksCont->GetArray()));
1065  if(!constituent)
1066  continue;
1067 
1068  // Plots on MC percentage in jet
1069  if(constituent->GetLabel() > 10000)
1070  {
1071  tracksFromMC++;
1072  ptFromMC += constituent->Pt();
1073  }
1074  tracksTotal++;
1075  ptTotal += constituent->Pt();
1076  }
1077  }
1078 
1079 
1080  trackRatio = 0.;
1081  if(tracksTotal)
1082  trackRatio = tracksFromMC/((Double_t)tracksTotal);
1083 
1085  ptTotal = jet->Pt() - fJetsCont->GetRhoVal()*jet->Area();
1086 
1087  ptRatio = 0.;
1088  if(ptTotal)
1089  ptRatio = ptFromMC/ptTotal;
1090 
1091 }
1092 
1093 //________________________________________________________________________
1095 {
1096  Double_t leadingVetoJetPt = -999.;
1097  AliEmcalJet* leadingVetoJet = 0;
1098  // Search for the 'leading' overlapping jet in the veto sample
1099  if(fJetVetoArray && jet)
1100  {
1101  for(Int_t j=0; j<fJetVetoArray->GetEntries(); j++)
1102  {
1103  UInt_t dummy = 0;
1104  AliEmcalJet* vetoJet = static_cast<AliEmcalJet*>(fJetVetoArray->At(j));
1105  // Check if veto jet is accepted
1106  if(!fJetsCont->AcceptJet(vetoJet , dummy))
1107  continue;
1108 
1109  // Check matching distance
1110  Double_t vetoPt = vetoJet->Pt() - vetoJet->Area()*fJetsCont->GetRhoVal();
1111  Double_t deltaEta = (jet->Eta()-vetoJet->Eta());
1112  Double_t deltaPhi = TMath::Min(TMath::Abs(jet->Phi()-vetoJet->Phi()),TMath::TwoPi() - TMath::Abs(jet->Phi()-vetoJet->Phi()));
1113  Double_t deltaR = TMath::Sqrt(deltaEta*deltaEta + deltaPhi*deltaPhi);
1114 
1115  if ((vetoPt > leadingVetoJetPt) && (deltaR <= fJetEmbeddingMaxDistance))
1116  {
1117  leadingVetoJetPt = vetoPt;
1118  leadingVetoJet = vetoJet;
1119  }
1120  }
1121  }
1122 
1123  return leadingVetoJet;
1124 }
1125 
1126 //________________________________________________________________________
1128 {
1129  Double_t leadingVetoJetPt = -999.;
1130  AliEmcalJet* leadingVetoJet = 0;
1131  // Search for the 'leading' jet in the veto sample
1132  if(fJetVetoArray)
1133  {
1134  for(Int_t j=0; j<fJetVetoArray->GetEntries(); j++)
1135  {
1136  UInt_t dummy = 0;
1137  AliEmcalJet* vetoJet = static_cast<AliEmcalJet*>(fJetVetoArray->At(j));
1138  // Check if veto jet is accepted
1139  if(!fJetsCont->AcceptJet(vetoJet , dummy))
1140  continue;
1141 
1142  // Check matching distance
1143  Double_t vetoPt = vetoJet->Pt() - vetoJet->Area()*fJetsCont->GetRhoVal();
1144  if (vetoPt > leadingVetoJetPt)
1145  {
1146  leadingVetoJetPt = vetoPt;
1147  leadingVetoJet = vetoJet;
1148  }
1149  }
1150  }
1151  return leadingVetoJet;
1152 }
1153 
1154 //________________________________________________________________________
1156 {
1157  // This is to use a full cone in phi even at the edges of phi (2pi -> 0) (0 -> 2pi)
1158  Double_t trackPhi = 0.0;
1159  if (track->Phi() > (TMath::TwoPi() - (radius-phi)))
1160  trackPhi = track->Phi() - TMath::TwoPi();
1161  else if (track->Phi() < (phi+radius - TMath::TwoPi()))
1162  trackPhi = track->Phi() + TMath::TwoPi();
1163  else
1164  trackPhi = track->Phi();
1165 
1166  if ( TMath::Abs(trackPhi-phi)*TMath::Abs(trackPhi-phi) + TMath::Abs(track->Eta()-eta)*TMath::Abs(track->Eta()-eta) <= radius*radius)
1167  return kTRUE;
1168 
1169  return kFALSE;
1170 }
1171 
1172 
1173 //________________________________________________________________________
1175 {
1176  // Calculate leading + subleading jet
1178  if(fJetOutputMode==6)
1180  else if(fJetOutputMode==4 || fJetOutputMode==5)
1181  GetMatchingJets();
1182 }
1183 
1184 //________________________________________________________________________
1185 void AliAnalysisTaskChargedJetsHadronCF::GetLeadingJets(const char* opt, AliEmcalJet*& jetLeading, AliEmcalJet*& jetSubLeading)
1186 {
1187  // Customized from AliJetContainer::GetLeadingJet()
1188  // Get the leading+subleading jet; if opt contains "rho" the sorting is according to pt-A*rho
1189 
1190  TString option(opt);
1191  option.ToLower();
1192 
1193  jetLeading = 0;
1194  jetSubLeading = 0;
1195 
1196  fJetsCont->ResetCurrentID();
1197  Double_t tmpLeadingPt = 0;
1198  Double_t tmpSubleadingPt = 0;
1199 
1200  if (option.Contains("rho")) {
1201  while (AliEmcalJet* jet = fJetsCont->GetNextAcceptJet()) {
1202  if ( (jet->Pt()-jet->Area()*fJetsCont->GetRhoVal()) > tmpLeadingPt )
1203  {
1204  jetSubLeading = jetLeading;
1205  jetLeading = jet;
1206  tmpSubleadingPt = tmpLeadingPt;
1207  tmpLeadingPt = jet->Pt()-jet->Area()*fJetsCont->GetRhoVal();
1208  }
1209  else if ( (jet->Pt()-jet->Area()*fJetsCont->GetRhoVal()) > tmpSubleadingPt )
1210  {
1211  jetSubLeading = jet;
1212  tmpSubleadingPt = jet->Pt()-jet->Area()*fJetsCont->GetRhoVal();
1213  }
1214  }
1215  }
1216  else {
1217  while (AliEmcalJet* jet = fJetsCont->GetNextAcceptJet()) {
1218  if ( (jet->Pt()) > tmpLeadingPt )
1219  {
1220  jetSubLeading = jetLeading;
1221  jetLeading = jet;
1222  tmpSubleadingPt = tmpLeadingPt;
1223  tmpLeadingPt = jet->Pt();
1224  }
1225  else if ( (jet->Pt()) > tmpSubleadingPt )
1226  {
1227  jetSubLeading = jet;
1228  tmpSubleadingPt = jet->Pt();
1229  }
1230  }
1231  }
1232 }
1233 
1234 //________________________________________________________________________
1235 void AliAnalysisTaskChargedJetsHadronCF::GetLeadingJetsInArray(TClonesArray* arr, const char* opt, AliEmcalJet*& jetLeading, AliEmcalJet*& jetSubLeading)
1236 {
1237  // Customized from AliJetContainer::GetLeadingJet()
1238  // Get the leading+subleading jet; if opt contains "rho" the sorting is according to pt-A*rho
1239 
1240  TString option(opt);
1241  option.ToLower();
1242 
1243  jetLeading = 0;
1244  jetSubLeading = 0;
1245  Double_t tmpLeadingPt = -999.;
1246  Double_t tmpSubleadingPt = -999.;
1247 
1248  for(Int_t i=0; i<arr->GetEntries(); i++)
1249  {
1250  AliEmcalJet* jet = static_cast<AliEmcalJet*>(arr->At(i));
1251  UInt_t dummy = 0;
1252  if(!fJetsCont->AcceptJet(jet , dummy))
1253  continue;
1254 
1255  Double_t jetPt = jet->Pt();
1256  if (option.Contains("rho"))
1257  jetPt -= jet->Area()*fJetsCont->GetRhoVal();
1258 
1259  if ( jetPt > tmpLeadingPt )
1260  {
1261  jetSubLeading = jetLeading;
1262  jetLeading = jet;
1263  tmpSubleadingPt = tmpLeadingPt;
1264  tmpLeadingPt = jetPt;
1265  }
1266  else if ( jetPt > tmpSubleadingPt )
1267  {
1268  jetSubLeading = jet;
1269  tmpSubleadingPt = jetPt;
1270  }
1271 
1272  }
1273 }
1274 
1275 
1276 //________________________________________________________________________
1278 {
1279  // Method for the correct logarithmic binning of histograms
1280  TAxis *axis = h->GetAxis(axisNumber);
1281  int bins = axis->GetNbins();
1282 
1283  Double_t from = axis->GetXmin();
1284  Double_t to = axis->GetXmax();
1285  Double_t *newBins = new Double_t[bins + 1];
1286 
1287  newBins[0] = from;
1288  Double_t factor = pow(to/from, 1./bins);
1289 
1290  for (int i = 1; i <= bins; i++) {
1291  newBins[i] = factor * newBins[i-1];
1292  }
1293  axis->Set(bins, newBins);
1294  delete [] newBins;
1295 }
1296 
1297 //________________________________________________________________________
1299 {
1300  std::vector<AliEmcalJet*>::iterator matchedJetFindResult = std::find(fMatchedJets.begin(), fMatchedJets.end(), jet);
1301  if(matchedJetFindResult == fMatchedJets.end())
1302  {
1303  AliError("Checked for a reference jet but it was not found. Check code.");
1304  return 0;
1305  }
1306 
1307  Int_t matchedJetIndex = (matchedJetFindResult - fMatchedJets.begin());
1308  AliEmcalJet* refJet = fMatchedJetsReference[matchedJetIndex];
1309  return refJet;
1310 }
1311 
1312 //________________________________________________________________________
1314 {
1315  TH1* tmpHist = static_cast<TH1*>(fOutput->FindObject(key));
1316  if(!tmpHist)
1317  {
1318  AliError(Form("Cannot find histogram <%s> ",key)) ;
1319  return;
1320  }
1321 
1322  tmpHist->Fill(x);
1323 }
1324 
1325 //________________________________________________________________________
1327 {
1328  TH1* tmpHist = static_cast<TH1*>(fOutput->FindObject(key));
1329  if(!tmpHist)
1330  {
1331  AliError(Form("Cannot find histogram <%s> ",key));
1332  return;
1333  }
1334 
1335  if (tmpHist->IsA()->GetBaseClass("TH1"))
1336  static_cast<TH1*>(tmpHist)->Fill(x,y); // Fill x with y
1337  else if (tmpHist->IsA()->GetBaseClass("TH2"))
1338  static_cast<TH2*>(tmpHist)->Fill(x,y); // Fill x,y with 1
1339 }
1340 
1341 //________________________________________________________________________
1343 {
1344  TH2* tmpHist = static_cast<TH2*>(fOutput->FindObject(key));
1345  if(!tmpHist)
1346  {
1347  AliError(Form("Cannot find histogram <%s> ",key));
1348  return;
1349  }
1350 
1351  tmpHist->Fill(x,y,add);
1352 }
1353 
1354 //________________________________________________________________________
1356 {
1357  TH3* tmpHist = static_cast<TH3*>(fOutput->FindObject(key));
1358  if(!tmpHist)
1359  {
1360  AliError(Form("Cannot find histogram <%s> ",key));
1361  return;
1362  }
1363 
1364  if(add)
1365  tmpHist->Fill(x,y,z,add);
1366  else
1367  tmpHist->Fill(x,y,z);
1368 }
1369 
1370 
1371 //________________________________________________________________________
1372 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)
1373 {
1374  T* tmpHist = new T(name, title, xBins, xMin, xMax);
1375 
1376  tmpHist->GetXaxis()->SetTitle(xTitle);
1377  tmpHist->GetYaxis()->SetTitle(yTitle);
1378  tmpHist->SetOption(options);
1379  tmpHist->SetMarkerStyle(kFullCircle);
1380  tmpHist->Sumw2();
1381 
1382  fOutput->Add(tmpHist);
1383 
1384  return tmpHist;
1385 }
1386 
1387 //________________________________________________________________________
1388 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)
1389 {
1390  T* tmpHist = new T(name, title, xBins, xMin, xMax, yBins, yMin, yMax);
1391  tmpHist->GetXaxis()->SetTitle(xTitle);
1392  tmpHist->GetYaxis()->SetTitle(yTitle);
1393  tmpHist->GetZaxis()->SetTitle(zTitle);
1394  tmpHist->SetOption(options);
1395  tmpHist->SetMarkerStyle(kFullCircle);
1396  tmpHist->Sumw2();
1397 
1398  fOutput->Add(tmpHist);
1399 
1400  return tmpHist;
1401 }
1402 
1403 //________________________________________________________________________
1404 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)
1405 {
1406  T* tmpHist = new T(name, title, xBins, xMin, xMax, yBins, yMin, yMax, zBins, zMin, zMax);
1407  tmpHist->GetXaxis()->SetTitle(xTitle);
1408  tmpHist->GetYaxis()->SetTitle(yTitle);
1409  tmpHist->GetZaxis()->SetTitle(zTitle);
1410  tmpHist->SetOption(options);
1411  tmpHist->SetMarkerStyle(kFullCircle);
1412  tmpHist->Sumw2();
1413 
1414  fOutput->Add(tmpHist);
1415 
1416  return tmpHist;
1417 }
1418 
1419 //________________________________________________________________________
1421 {
1422  // Called once at the end of the analysis.
1423 }
1424 
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:110
void AddJetConstituent(Float_t eta, Float_t phi, Float_t pt, Short_t charge, Short_t pid=0, Float_t vx=0, Float_t vy=0, Float_t vz=0)
Double_t Area() const
Definition: AliEmcalJet.h:117
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:26
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:108
long long Long64_t
Definition: External.C:43
Double_t Phi() const
Definition: AliEmcalJet.h:104
void * fJetsTreeBuffer
! buffer for one jet (that will be saved to the tree)
Container with name, TClonesArray and cuts for particles.
Support task for (charged) jet-hadron correlations.
Int_t fPythiaExtractionMode
Mode which PYTHIA-jets to extract for fJetOutputMode==6: 0: all, 1: quark-jets, 2: gluon jets...
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
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)
Int_t GetPartonFlag7() const
AliEmcalJet * fSubleadingJet
! subleading jet (calculated event-by-event)
Simple class containing basic information for a jet.
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:147
Container class of cuts for AliAnalysisTaskChargedJetsHadronCF.
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:126
UShort_t T(UShort_t m, UShort_t t)
Definition: RingBits.C:60
Bool_t fPythiaExtractionUseHadronMatching
Use hadron matching to select/label jets (works for b,c,s and others)
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
Int_t GetPartonFlag6() const
void GetTrackMCRatios(AliEmcalJet *jet, AliEmcalJet *mcJet, Double_t &trackRatio, Double_t &ptRatio)
AliParticleContainer * GetParticleContainer() const
Bool_t fJetEmbeddingCreatePtPlotPerCut
create TH3 per cut or only once
int Int_t
Definition: External.C:63
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.
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.
Double_t fExtractionPercentage
percentage that is recorded
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)
Double_t fExtractionMinPt
minimum pt of recorded jets
AliEmcalJet * GetNextAcceptJet()
void AddJetToOutputArray(AliEmcalJet *jet, Int_t arrayIndex, Int_t &jetsAlreadyInArray)
TTree * fJetsTree
! Jets that will be saved to a tree (optionally)
Double_t Pt() const
Definition: AliEmcalJet.h:96
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.
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
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
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:44
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
Double_t fExtractionMaxPt
maximum pt of recorded jets
virtual AliVParticle * GetNextAcceptParticle()
bool Bool_t
Definition: External.C:53
void CalculateJetType(AliEmcalJet *jet, Int_t &typeIC, Int_t &typeHM)
Double_t fJetEmbeddingMaxDistance
Max distance allowed to accept an embedded jet.
void FillHistogramsJets(AliEmcalJet *jet, const char *cutName)
Double_t fHadronMatchingRadius
Radius used in the hadron matching (MC jet extraction)
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 fConstPtFilterBit
For const pt plot, filter bit.