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