AliPhysics  7273240 (7273240)
 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 
30 #include "AliVTrack.h"
31 #include "AliVHeader.h"
32 #include "AliEmcalJet.h"
33 #include "AliRhoParameter.h"
34 #include "AliLog.h"
35 #include "AliJetContainer.h"
36 #include "AliTrackContainer.h"
37 #include "AliAODTrack.h"
38 #include "AliPicoTrack.h"
39 #include "AliVParticle.h"
40 #include "TRandom3.h"
42 
44 
48 //________________________________________________________________________
50 {
51 // dummy destructor
52 }
53 
57 //________________________________________________________________________
59 {
60 // dummy destructor
61 }
62 
66 //________________________________________________________________________
69  fJetsCont(0),
70  fTracksCont(0),
71  fJetsTree(0),
72  fJetsTreeBuffer(0),
73  fExtractionPercentage(0),
74  fExtractionMinPt(0),
75  fExtractionMaxPt(0),
76  fEventExtractionPercentage(0),
77  fEventExtractionMinJetPt(0),
78  fEventExtractionMaxJetPt(0),
79  fNumberOfCentralityBins(10),
80  fJetsOutput(),
81  fTracksOutput(),
82  fJetParticleArrayName("JetsDPhiBasicParticles"),
83  fTrackParticleArrayName(""),
84  fJetMatchingArray(),
85  fJetMatchingArrayName(""),
86  fJetMatchingMaxDistance(0.3),
87  fJetMatchingMinSharedFraction(0.5),
88  fJetMatchingMaxSharedFraction(1.5),
89  fJetMatchingMaxEmbeddingOffset(20.),
90  fJetMatchingMinPt(0.0),
91  fJetMatchingMaxPt(999.0),
92  fJetMatchingUseOnlyNLeading(0),
93  fJetVetoArray(),
94  fJetVetoArrayName(""),
95  fJetVetoMinPt(0),
96  fJetVetoMaxPt(0),
97  fMatchedJets(),
98  fRandom(0),
99  fJetOutputMode(0),
100  fMinFakeFactorPercentage(0),
101  fMaxFakeFactorPercentage(0),
102  fPythiaExtractionMode(0),
103  fEventCriteriumMode(0),
104  fEventCriteriumMinBackground(0),
105  fEventCriteriumMaxBackground(0),
106  fEventCriteriumMinLeadingJetPt(0),
107  fEventCriteriumMinSubleadingJetPt(0),
108  fEventCriteriumMinJetDeltaPhi(0),
109  fLeadingJet(0),
110  fSubleadingJet(0),
111  fInitialPartonMatchedJet1(0),
112  fInitialPartonMatchedJet2(0),
113  fAcceptedJets(0),
114  fAcceptedTracks(0)
115 {
116  // Default constructor.
117  SetMakeGeneralHistograms(kTRUE);
118  fRandom = new TRandom3(0);
119 }
120 
121 
122 //________________________________________________________________________
124  AliAnalysisTaskEmcalJet(name, kTRUE),
125  fJetsCont(0),
126  fTracksCont(0),
127  fJetsTree(0),
128  fJetsTreeBuffer(0),
129  fExtractionPercentage(0),
130  fExtractionMinPt(0),
131  fExtractionMaxPt(0),
132  fEventExtractionPercentage(0),
133  fEventExtractionMinJetPt(0),
134  fEventExtractionMaxJetPt(0),
135  fNumberOfCentralityBins(10),
136  fJetsOutput(),
137  fTracksOutput(),
138  fJetParticleArrayName("JetsDPhiBasicParticles"),
139  fTrackParticleArrayName(""),
140  fJetMatchingArray(),
141  fJetMatchingArrayName(""),
142  fJetMatchingMaxDistance(0.3),
143  fJetMatchingMinSharedFraction(0.5),
144  fJetMatchingMaxSharedFraction(1.5),
145  fJetMatchingMaxEmbeddingOffset(20.),
146  fJetMatchingMinPt(0.0),
147  fJetMatchingMaxPt(999.0),
148  fJetMatchingUseOnlyNLeading(0),
149  fJetVetoArray(),
150  fJetVetoArrayName(""),
151  fJetVetoMinPt(0),
152  fJetVetoMaxPt(0),
153  fMatchedJets(),
154  fRandom(0),
155  fJetOutputMode(0),
156  fMinFakeFactorPercentage(0),
157  fMaxFakeFactorPercentage(0),
158  fPythiaExtractionMode(0),
159  fEventCriteriumMode(0),
160  fEventCriteriumMinBackground(0),
161  fEventCriteriumMaxBackground(0),
162  fEventCriteriumMinLeadingJetPt(0),
163  fEventCriteriumMinSubleadingJetPt(0),
164  fEventCriteriumMinJetDeltaPhi(0),
165  fLeadingJet(0),
166  fSubleadingJet(0),
167  fInitialPartonMatchedJet1(0),
168  fInitialPartonMatchedJet2(0),
169  fAcceptedJets(0),
170  fAcceptedTracks(0)
171 {
172  // Constructor
174  fRandom = new TRandom3(0);
175 }
176 
177 //________________________________________________________________________
179 {
180  // Destructor.
181 }
182 
183 //________________________________________________________________________
185 {
187 
188  // ### Basic container settings
190  if(fJetsCont) { //get particles connected to jets
191  fJetsCont->PrintCuts();
193  } else { //no jets, just analysis tracks
195  }
196  if(fTracksCont) fTracksCont->SetClassName("AliVTrack");
197 
198  // ### Create all histograms
199 
200  // Change the event rejection histogram -> Add a custom value
201  fHistEventRejection->GetXaxis()->SetBinLabel(14,"JetCrit");
202 
203  // Track QA plots
204  AddHistogram2D<TH2D>("hTrackPt", "Tracks p_{T} distribution", "", 300, 0., 300., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)", "Centrality", "dN^{Tracks}/dp_{T}");
205  AddHistogram2D<TH2D>("hTrackPhi", "Track angular distribution in #phi", "LEGO2", 180, 0., 2*TMath::Pi(), fNumberOfCentralityBins, 0, 100, "#phi", "Centrality", "dN^{Tracks}/(d#phi)");
206  AddHistogram2D<TH2D>("hTrackEta", "Track angular distribution in #eta", "LEGO2", 100, -2.5, 2.5, fNumberOfCentralityBins, 0, 100, "#eta", "Centrality", "dN^{Tracks}/(d#eta)");
207  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");
208 
209  AddHistogram2D<TH2D>("hLeadingTrackPt", "Leading tracks p_{T} distribution", "", 300, 0., 300., fNumberOfCentralityBins, 0, 100, "p_{T} (GeV/c)", "Centrality", "dN^{Tracks}/dp_{T}");
210  AddHistogram2D<TH2D>("hLeadingTrackPhi", "Leading tracks angular distribution in #phi", "LEGO2", 180, 0., 2*TMath::Pi(), fNumberOfCentralityBins, 0, 100, "#phi", "Centrality", "dN^{Tracks}/(d#phi)");
211  AddHistogram2D<TH2D>("hLeadingTrackEta", "Leading tracks angular distribution in #eta", "LEGO2", 100, -2.5, 2.5, fNumberOfCentralityBins, 0, 100, "#eta", "Centrality", "dN^{Tracks}/(d#eta)");
212  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");
213 
214  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})");
215  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})");
216 
217 
218  // Jet QA plots
219  AddHistogram2D<TH2D>("hJetPtRaw", "Jets p_{T} distribution (no bgrd. corr.)", "", 300, 0., 300., fNumberOfCentralityBins, 0, 100, "p_{T, jet} (GeV/c)", "Centrality", "dN^{Jets}/dp_{T}");
220  AddHistogram2D<TH2D>("hJetPt", "Jets p_{T} distribution (background subtracted)", "", 400, -100., 300., fNumberOfCentralityBins, 0, 100, "p_{T, jet} (GeV/c)", "Centrality", "dN^{Jets}/dp_{T}");
221  AddHistogram2D<TH2D>("hJetPhi", "Jet angular distribution #phi", "LEGO2", 180, 0., 2*TMath::Pi(), fNumberOfCentralityBins, 0, 100, "#phi", "Centrality", "dN^{Jets}/d#phi");
222  AddHistogram2D<TH2D>("hJetEta", "Jet angular distribution #eta", "LEGO2", 100, -2.5, 2.5, fNumberOfCentralityBins, 0, 100, "#eta","Centrality","dN^{Jets}/d#eta");
223  AddHistogram2D<TH2D>("hJetPhiPt", "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}");
224  AddHistogram2D<TH2D>("hJetEtaPt", "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}");
225  AddHistogram2D<TH2D>("hJetPhiEta", "Jet angular distribution #phi/#eta", "COLZ", 180, 0., 2*TMath::Pi(), 100, -2.5, 2.5, "#phi", "#eta", "dN^{Jets}/d#phi d#eta");
226  AddHistogram2D<TH2D>("hJetArea", "Jet area", "LEGO2", 200, 0., 2., fNumberOfCentralityBins, 0, 100, "Jet A", "Centrality", "dN^{Jets}/dA");
227  AddHistogram2D<TH2D>("hJetAreaPt", "Jet area vs. p_{T}", "LEGO2", 200, 0., 2., 400, -100., 300., "Jet A", "p_{T, jet} (GeV/c)", "dN^{Jets}/dA dp_{T}");
228  AddHistogram2D<TH2D>("hJetPtLeadingHadron", "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}");
229 
230  AddHistogram2D<TH2D>("hJetConstituentPt_Cent0_100", "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}");
231  AddHistogram2D<TH2D>("hJetConstituentPt_Cent0_10", "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}");
232 
233  AddHistogram2D<TH2D>("hJetConstituentCount_Cent0_100", "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}");
234  AddHistogram2D<TH2D>("hJetConstituentCount_Cent0_10", "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}");
235 
236  // Embedding plots
237  if(fJetOutputMode == 4 || fJetOutputMode == 5)
238  {
239  AddHistogram2D<TH2D>("hEmbeddingDeltaR", "Matched jet #Delta R distribution", "", 200, -50., 150., 100, 0, 1.0, "p_{T, jet} (GeV/c)", "#Delta R", "dN^{Matched}/dp_{T}dR");
240  AddHistogram2D<TH2D>("hEmbeddingDeltaEta", "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");
241  AddHistogram2D<TH2D>("hEmbeddingDeltaPhi", "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");
242  AddHistogram2D<TH2D>("hEmbeddingPtCorr", "Matched jet p_{T} distributions", "", 200, -50., 150., 200, -50., 150., "p_{T, MC jet} (GeV/c)", "p_{T, emb} (GeV/c)", "dN^{Matched}/dp_{T}d#Delta p_{T}");
243 
244  AddHistogram1D<TH1D>("hEmbeddingJetPt", "Embedded jets p_{T} distribution", "", 200, -50., 150., "p_{T, jet} (GeV/c)", "dN/dp_{T}");
245  AddHistogram2D<TH2D>("hEmbeddingJetPhiEta", "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");
246  }
247 
248  // Random cone plots
249  AddHistogram2D<TH2D>("hRandomConePt", "Random cone p_{T} distribution", "", 400, -100., 300., fNumberOfCentralityBins, 0, 100, "p_{T, cone} (GeV/c)", "Centrality", "dN^{Tracks}/dp_{T}");
250  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}");
251  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}");
252  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}");
253 
254  // Leading/subleading, background ...
255 
256  AddHistogram2D<TH2D>("hLeadingJetPtRaw", "Jets p_{T} distribution (no bgrd. corr.)", "", 300, 0., 300., fNumberOfCentralityBins, 0, 100, "p_{T, jet} (GeV/c)", "Centrality", "dN^{Jets}/dp_{T}");
257  AddHistogram2D<TH2D>("hLeadingJetPt", "Jets p_{T} distribution (background subtracted)", "", 400, -100., 300., fNumberOfCentralityBins, 0, 100, "p_{T, jet} (GeV/c)", "Centrality", "dN^{Jets}/dp_{T}");
258  AddHistogram2D<TH2D>("hLeadingJetPhi", "Jet angular distribution #phi", "LEGO2", 180, 0., 2*TMath::Pi(), fNumberOfCentralityBins, 0, 100, "#phi", "Centrality", "dN^{Jets}/d#phi");
259  AddHistogram2D<TH2D>("hLeadingJetEta", "Jet angular distribution #eta", "LEGO2", 100, -2.5, 2.5, fNumberOfCentralityBins, 0, 100, "#eta","Centrality","dN^{Jets}/d#eta");
260  AddHistogram2D<TH2D>("hLeadingJetPhiPt", "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}");
261  AddHistogram2D<TH2D>("hLeadingJetEtaPt", "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}");
262  AddHistogram2D<TH2D>("hLeadingJetPhiEta", "Jet angular distribution #phi/#eta", "COLZ", 180, 0., 2*TMath::Pi(), 100, -2.5, 2.5, "#phi", "#eta", "dN^{Jets}/d#phi d#eta");
263  AddHistogram2D<TH2D>("hLeadingJetArea", "Jet area", "LEGO2", 200, 0., 2., fNumberOfCentralityBins, 0, 100, "Jet A", "Centrality", "dN^{Jets}/dA");
264  AddHistogram2D<TH2D>("hLeadingJetAreaPt", "Jet area vs. p_{T}", "LEGO2", 200, 0., 2., 400, -100., 300., "Jet A", "p_{T, jet} (GeV/c)", "dN^{Jets}/dA dp_{T}");
265  AddHistogram2D<TH2D>("hLeadingJetPtLeadingHadron", "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}");
266 
267  AddHistogram2D<TH2D>("hSubleadingJetPtRaw", "Jets p_{T} distribution (no bgrd. corr.)", "", 300, 0., 300., fNumberOfCentralityBins, 0, 100, "p_{T, jet} (GeV/c)", "Centrality", "dN^{Jets}/dp_{T}");
268  AddHistogram2D<TH2D>("hSubleadingJetPt", "Jets p_{T} distribution (background subtracted)", "", 400, -100., 300., fNumberOfCentralityBins, 0, 100, "p_{T, jet} (GeV/c)", "Centrality", "dN^{Jets}/dp_{T}");
269  AddHistogram2D<TH2D>("hSubleadingJetPhi", "Jet angular distribution #phi", "LEGO2", 180, 0., 2*TMath::Pi(), fNumberOfCentralityBins, 0, 100, "#phi", "Centrality", "dN^{Jets}/d#phi");
270  AddHistogram2D<TH2D>("hSubleadingJetEta", "Jet angular distribution #eta", "LEGO2", 100, -2.5, 2.5, fNumberOfCentralityBins, 0, 100, "#eta","Centrality","dN^{Jets}/d#eta");
271  AddHistogram2D<TH2D>("hSubleadingJetPhiPt", "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}");
272  AddHistogram2D<TH2D>("hSubleadingJetEtaPt", "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}");
273  AddHistogram2D<TH2D>("hSubleadingJetPhiEta", "Jet angular distribution #phi/#eta", "COLZ", 180, 0., 2*TMath::Pi(), 100, -2.5, 2.5, "#phi", "#eta", "dN^{Jets}/d#phi d#eta");
274  AddHistogram2D<TH2D>("hSubleadingJetArea", "Jet area", "LEGO2", 200, 0., 2., fNumberOfCentralityBins, 0, 100, "Jet A", "Centrality", "dN^{Jets}/dA");
275  AddHistogram2D<TH2D>("hSubleadingJetAreaPt", "Jet area vs. p_{T}", "LEGO2", 200, 0., 2., 400, -100., 300., "Jet A", "p_{T, jet} (GeV/c)", "dN^{Jets}/dA dp_{T}");
276  AddHistogram2D<TH2D>("hSubleadingJetPtLeadingHadron", "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}");
277 
278  AddHistogram2D<TH2D>("hTrackCount", "Number of tracks in acceptance vs. centrality", "LEGO2", 500, 0., 5000., fNumberOfCentralityBins, 0, 100, "N tracks","Centrality", "dN^{Events}/dN^{Tracks}");
279  AddHistogram2D<TH2D>("hJetCount", "Number of jets in acceptance vs. centrality", "LEGO2", 100, 0., 100., fNumberOfCentralityBins, 0, 100, "N Jets","Centrality", "dN^{Events}/dN^{Jets}");
280  AddHistogram2D<TH2D>("hFakeFactor", "Fake factor distribution", "LEGO2", 1000, 0., 100., fNumberOfCentralityBins, 0, 100, "Fake factor","Centrality", "dN^{Jets}/df");
281  AddHistogram2D<TH2D>("hFakeFactorJetPt_Cent0_100", "Fake factor distribution vs. jet p_{T}", "LEGO2", 1000, 0., 100., 400, -100., 300., "Fake factor","Jet p_{T} (GeV/c)", "dN^{Jets}/df");
282  AddHistogram2D<TH2D>("hFakeFactorJetPt_Cent0_10", "Fake factor distribution vs. jet p_{T}", "LEGO2", 1000, 0., 100., 400, -100., 300., "Fake factor","Jet p_{T} (GeV/c)", "dN^{Jets}/df");
283 
284  AddHistogram2D<TH2D>("hBackgroundPt", "Background p_{T} distribution", "", 150, 0., 150., fNumberOfCentralityBins, 0, 100, "Background p_{T} (GeV/c)", "Centrality", "dN^{Events}/dp_{T}");
285  AddHistogram2D<TH2D>("hBackgroundPtJetPt_Cent0_100", "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}");
286  AddHistogram2D<TH2D>("hBackgroundPtJetPt_Cent0_10", "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}");
287  AddHistogram2D<TH2D>("hBackgroundPtConstCount_Cent0_100", "Background p_{T} distribution vs. const. count", "", 150, 0., 150., 200, 0., 200., "Background p_{T} (GeV/c)", "Count", "dN^{Events}/dp_{T}");
288  AddHistogram2D<TH2D>("hBackgroundPtConstCount_Cent0_10", "Background p_{T} distribution vs. const. count", "", 150, 0., 150., 200, 0., 200., "Background p_{T} (GeV/c)", "Count", "dN^{Events}/dp_{T}");
289 
290  PostData(1, fOutput); // Post data for ALL output slots > 0 here.
291 }
292 
293 
294 //________________________________________________________________________
296 
298 
299  // ### Add the jets as basic correlation particles to the event
300  if (!(fInputEvent->FindListObject(Form("%s", fJetParticleArrayName.Data()))))
301  {
302  fJetsOutput = new TClonesArray("AliPicoTrack");
303  fJetsOutput->SetName(fJetParticleArrayName.Data());
304  fInputEvent->AddObject(fJetsOutput);
305  }
306  else
307  AliError(Form("%s: Object with name %s already in event!", GetName(), Form("%s", fJetParticleArrayName.Data())));
308 
309  // ### Add the tracks as basic correlation particles to the event (optional)
310  if(fTrackParticleArrayName != "")
311  {
312  if (!(fInputEvent->FindListObject(Form("%s", fTrackParticleArrayName.Data()))))
313  {
314  fTracksOutput = new TClonesArray("AliPicoTrack");
315  fTracksOutput->SetName(fTrackParticleArrayName.Data());
316  fInputEvent->AddObject(fTracksOutput);
317  }
318  else
319  AliError(Form("%s: Object with name %s already in event!", GetName(), Form("%s", fTrackParticleArrayName.Data())));
320  }
321 
322  // ### Import generated jets from toymodel for matching (optional)
323  if(fJetMatchingArrayName != "")
324  {
325  fJetMatchingArray = static_cast<TClonesArray*>(InputEvent()->FindListObject(Form("%s", fJetMatchingArrayName.Data())));
326  if(!fJetMatchingArray)
327  AliFatal(Form("Importing jets for matching failed! Array '%s' not found!", fJetMatchingArrayName.Data()));
328  }
329  else if(fJetOutputMode==4 || fJetOutputMode==5)
330  AliFatal(Form("fJetMatchingArrayName must be set in jet output mode 4 or 5."));
331 
332  // ### Import veto jets for matching (optional)
333  if(fJetVetoArrayName != "")
334  {
335  fJetVetoArray = static_cast<TClonesArray*>(InputEvent()->FindListObject(Form("%s", fJetVetoArrayName.Data())));
336  if(!fJetVetoArray)
337  AliFatal(Form("Importing jets for veto failed! Array '%s' not found!", fJetVetoArrayName.Data()));
338  }
339 
340  // ### Jets tree (optional)
342  {
343  fJetsTree = new TTree("ExtractedJets", "ExtractedJets");
344  fJetsTree->Branch("Jets", "AliBasicJet", &fJetsTreeBuffer, 1000);
345  fOutput->Add(fJetsTree);
346  }
347 
348 
349 }
350 
351 //________________________________________________________________________
353 {
354 
355  // In case of special selection criteria, trigger on certain events
356  if(fEventCriteriumMode==0) // "minimum bias"
357  {
358  // do nothing
359  }
360  else if(fEventCriteriumMode==1) // background constraints
361  {
363  {
364  fHistEventRejection->Fill("JetCrit", 1);
365  return kFALSE;
366  }
367  }
368  else if(fEventCriteriumMode==2) // Minimum leading jet pT
369  {
370  if(fLeadingJet)
371  {
373  {
374  fHistEventRejection->Fill("JetCrit", 1);
375  return kFALSE;
376  }
377  }
378  }
379  else if(fEventCriteriumMode==3) // Simple dijet trigger
380  {
382  {
385  {
386  fHistEventRejection->Fill("JetCrit", 1);
387  return kFALSE;
388  }
389  else // dijet pT fulfilled, check back-to-back criterium
390  {
391  Double_t deltaPhi = TMath::Min(TMath::Abs(fLeadingJet->Phi()-fSubleadingJet->Phi()),TMath::TwoPi() - TMath::Abs(fLeadingJet->Phi()-fSubleadingJet->Phi()));
392  if(deltaPhi <= fEventCriteriumMinJetDeltaPhi)
393  {
394  fHistEventRejection->Fill("JetCrit", 1);
395  return kFALSE;
396  }
397  }
398  }
399  }
400  return kTRUE;
401 }
402 
403 //________________________________________________________________________
405 {
406  if(fJetOutputMode==3) // output leading&subleading jet
407  {
408  if((jet!=fLeadingJet) && (jet!=fSubleadingJet))
409  return kFALSE;
410  }
411  else if(fJetOutputMode==1) // output the leading jet
412  {
413  if(jet!=fLeadingJet)
414  return kFALSE;
415  }
416  else if(fJetOutputMode==2) // output the subleading jet
417  {
418  if(jet!=fSubleadingJet)
419  return kFALSE;
420  }
421  else if(fJetOutputMode==6)
422  {
424  return kFALSE;
425  }
426 
427  if(fJetOutputMode==4) // matching jets only
428  return (std::find(fMatchedJets.begin(), fMatchedJets.end(), jet) != fMatchedJets.end());
429  else if(fJetOutputMode==5) // non-matching jets only
430  return (std::find(fMatchedJets.begin(), fMatchedJets.end(), jet) == fMatchedJets.end());
431 
432  return kTRUE;
433 }
434 
435 //________________________________________________________________________
437 {
438  // All jets
439  FillHistogram("hJetPtRaw", jet->Pt(), fCent);
440  FillHistogram("hJetPt", jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), fCent);
441  FillHistogram("hJetPhi", jet->Phi(), fCent);
442  FillHistogram("hJetEta", jet->Eta(), fCent);
443  FillHistogram("hJetEtaPt", jet->Eta(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
444  FillHistogram("hJetPhiPt", jet->Phi(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
445  FillHistogram("hJetPhiEta", jet->Phi(), jet->Eta());
446  FillHistogram("hJetArea", jet->Area(), fCent);
447  FillHistogram("hJetAreaPt", jet->Area(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
448  FillHistogram("hJetPtLeadingHadron", jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), fJetsCont->GetLeadingHadronPt(jet));
449 
450  FillHistogram("hBackgroundPtJetPt_Cent0_100", fJetsCont->GetRhoVal(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
451  if( (fCent >= 0) && (fCent < 10) )
452  FillHistogram("hBackgroundPtJetPt_Cent0_10", fJetsCont->GetRhoVal(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
453 
454  // Fake jet rejection (0810.1219)
455  Double_t fakeFactor = CalculateFakeFactor(jet);
456  FillHistogram("hFakeFactor", fakeFactor, fCent);
457  FillHistogram("hFakeFactorJetPt_Cent0_100", fakeFactor, jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
458  if( (fCent >= 0) && (fCent < 10) )
459  FillHistogram("hFakeFactorJetPt_Cent0_10", fakeFactor, jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
460 
461 
462  // Leading jet plots
463  if(jet==fLeadingJet)
464  {
465  FillHistogram("hLeadingJetPtRaw", jet->Pt(), fCent);
466  FillHistogram("hLeadingJetPt", jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), fCent);
467  FillHistogram("hLeadingJetPhi", jet->Phi(), fCent);
468  FillHistogram("hLeadingJetEta", jet->Eta(), fCent);
469  FillHistogram("hLeadingJetEtaPt", jet->Eta(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
470  FillHistogram("hLeadingJetPhiPt", jet->Phi(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
471  FillHistogram("hLeadingJetPhiEta", jet->Phi(), jet->Eta());
472  FillHistogram("hLeadingJetArea", jet->Area(), fCent);
473  FillHistogram("hLeadingJetAreaPt", jet->Area(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
474  FillHistogram("hLeadingJetPtLeadingHadron", jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), fJetsCont->GetLeadingHadronPt(jet));
475  }
476 
477  // Subleading jet plot
478  else if(jet==fSubleadingJet)
479  {
480  FillHistogram("hSubleadingJetPtRaw", jet->Pt(), fCent);
481  FillHistogram("hSubleadingJetPt", jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), fCent);
482  FillHistogram("hSubleadingJetPhi", jet->Phi(), fCent);
483  FillHistogram("hSubleadingJetEta", jet->Eta(), fCent);
484  FillHistogram("hSubleadingJetEtaPt", jet->Eta(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
485  FillHistogram("hSubleadingJetPhiPt", jet->Phi(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
486  FillHistogram("hSubleadingJetPhiEta", jet->Phi(), jet->Eta());
487  FillHistogram("hSubleadingJetArea", jet->Area(), fCent);
488  FillHistogram("hSubleadingJetAreaPt", jet->Area(), jet->Pt() - fJetsCont->GetRhoVal()*jet->Area());
489  FillHistogram("hSubleadingJetPtLeadingHadron", jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), fJetsCont->GetLeadingHadronPt(jet));
490  }
491 }
492 
493 //________________________________________________________________________
495 {
496  FillHistogram("hTrackPt", track->Pt(), fCent);
497  FillHistogram("hTrackPhi", track->Phi(), fCent);
498  FillHistogram("hTrackEta", track->Eta(), fCent);
499  FillHistogram("hTrackEtaPt", track->Eta(), track->Pt());
500  FillHistogram("hTrackPhiPt", track->Phi(), track->Pt());
501  FillHistogram("hTrackPhiEta", track->Phi(), track->Eta());
502 }
503 
504 //________________________________________________________________________
506 {
507  // Loop over all jet constituents
508  for(Int_t i = 0; i < jet->GetNumberOfTracks(); i++)
509  {
510  AliVParticle* constituent = static_cast<AliVParticle*>(jet->TrackAt(i, fTracksCont->GetArray()));
511  if(!constituent)
512  continue;
513 
514  // Fill jet constituent plots
515  FillHistogram("hJetConstituentPt_Cent0_100", jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), constituent->Pt());
516  if( (fCent >= 0) && (fCent < 10) )
517  FillHistogram("hJetConstituentPt_Cent0_10", jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), constituent->Pt());
518 
519  }
520 
521  FillHistogram("hJetConstituentCount_Cent0_100", jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), jet->GetNumberOfTracks());
522  if( (fCent >= 0) && (fCent < 10) )
523  FillHistogram("hJetConstituentCount_Cent0_10", jet->Pt() - fJetsCont->GetRhoVal()*jet->Area(), jet->GetNumberOfTracks());
524 
525  FillHistogram("hBackgroundPtConstCount_Cent0_100", fJetsCont->GetRhoVal(), jet->GetNumberOfTracks());
526  if( (fCent >= 0) && (fCent < 10) )
527  FillHistogram("hBackgroundPtConstCount_Cent0_10", fJetsCont->GetRhoVal(), jet->GetNumberOfTracks());
528 }
529 
530 //________________________________________________________________________
532 {
533  Double_t tmpPt = jet->Pt() - fJetsCont->GetRhoVal()*jet->Area();
534  new ((*fJetsOutput)[fAcceptedJets]) AliPicoTrack(tmpPt, jet->Eta(), jet->Phi(), jet->Charge(), 0, 0);
535  fAcceptedJets++;
536 }
537 
538 //________________________________________________________________________
540 {
541  if(fTrackParticleArrayName != "")
542  {
543  new ((*fTracksOutput)[fAcceptedTracks]) AliPicoTrack(track->Pt(), track->Eta(), track->Phi(), track->Charge(), 0, 0); // only Pt,Eta,Phi are interesting for correlations;
544  fAcceptedTracks++;
545  }
546 }
547 
548 //________________________________________________________________________
550 {
551  // Check pT threshold
552  if( ((jet->Pt()-jet->Area()*fJetsCont->GetRhoVal()) < fExtractionMinPt) || ((jet->Pt()-jet->Area()*fJetsCont->GetRhoVal()) >= fExtractionMaxPt) )
553  return;
554 
555  // Discard jets statistically
556  if(fRandom->Rndm() >= fExtractionPercentage)
557  return;
558 
559  AliVHeader* eventIDHeader = InputEvent()->GetHeader();
560  Long64_t eventID = 0;
561  if(eventIDHeader)
562  eventID = eventIDHeader->GetEventIdAsLong();
563 
564  // if only the two initial collision partons will be added, get PYTHIA info on them
565  Int_t partid = 0;
566  if(fJetOutputMode==6)
567  {
568  if(!fPythiaInfo)
569  AliError("fPythiaInfo object not available. Is it activated with SetGeneratePythiaInfoObject()?");
570  else if(jet==fInitialPartonMatchedJet1)
571  partid = fPythiaInfo->GetPartonFlag6();
572  else if (jet==fInitialPartonMatchedJet2)
573  partid = fPythiaInfo->GetPartonFlag7();
574 
575  // If fPythiaExtractionMode is set, only extract certain jets
576  if( (fPythiaExtractionMode==1) && not (partid>=1 && partid<=6)) // all quark-jet extraction
577  return;
578  else if( (fPythiaExtractionMode==2) && not (partid==21)) // gluon-jet extraction
579  return;
580  else if( (fPythiaExtractionMode<0) && (fPythiaExtractionMode!=-partid) ) // custom type jet extraction by given a negative number
581  return;
582  }
583 
584  AliBasicJet basicJet(jet->Eta(), jet->Phi(), jet->Pt(), jet->Charge(), fJetsCont->GetJetRadius(), jet->Area(), partid, fJetsCont->GetRhoVal(), eventID, fCent);
585  // Add constituents
586  for(Int_t i = 0; i < jet->GetNumberOfTracks(); i++)
587  {
588  AliVParticle* particle = static_cast<AliVParticle*>(jet->TrackAt(i, fTracksCont->GetArray()));
589  if(!particle) continue;
590 
591  AliAODTrack* aodtrack = static_cast<AliAODTrack*>(jet->TrackAt(i, fTracksCont->GetArray()));
592  Int_t constid = 9; // 9 mean unknown
593  if(fJetOutputMode==6)
594  {
595  // Use same convention as PID in AODs
596  if(TMath::Abs(particle->PdgCode()) == 2212) // proton
597  constid = 4;
598  else if (TMath::Abs(particle->PdgCode()) == 211) // pion
599  constid = 2;
600  else if (TMath::Abs(particle->PdgCode()) == 321) // kaon
601  constid = 3;
602  else if (TMath::Abs(particle->PdgCode()) == 11) // electron
603  constid = 0;
604  else if (TMath::Abs(particle->PdgCode()) == 13) // muon
605  constid = 1;
606  }
607  else if (aodtrack)
608  constid = aodtrack->GetMostProbablePID();
609 
610  basicJet.AddJetConstituent(particle->Eta(), particle->Phi(), particle->Pt(), particle->Charge(), constid);
611  }
612  if(std::find(fMatchedJets.begin(), fMatchedJets.end(), jet) != fMatchedJets.end()) // set the true pT from the matched jets (only possible in modes 4 & 7)
613  basicJet.SetTruePt(fMatchedJetsReference[std::find(fMatchedJets.begin(), fMatchedJets.end(), jet)-fMatchedJets.begin()]->Pt());
614 
615  fJetsTreeBuffer = &basicJet;
616  fJetsTree->Fill();
617 }
618 
619 //________________________________________________________________________
621 {
622  // Check jet pT threshold
624  return;
625 
626  // Discard jets statistically
627  if(fRandom->Rndm() >= fEventExtractionPercentage)
628  return;
629 
630  static Int_t numSavedEvents = 0;
631  numSavedEvents++;
632 
633 
634  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");
635  fTracksCont->ResetCurrentID();
636  while(AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle()))
637  FillHistogram(Form("Event%i", numSavedEvents), track->Phi(), track->Eta(), track->Pt());
638 
639 }
640 
641 //________________________________________________________________________
643 {
645 
647  return kFALSE;
648 
649  // ####### Jet loop
650  fAcceptedJets = 0;
651  fJetsCont->ResetCurrentID();
652  while(AliEmcalJet *jet = fJetsCont->GetNextAcceptJet())
653  {
654  if(!IsJetSelected(jet))
655  continue;
656 
657  // Jet plots
658  FillHistogramsJets(jet);
660 
661  // Add jet to output array
663  AddJetToTree(jet);
664  AddJetToOutputArray(jet);
665  }
666 
667  // ####### Particle loop
668  // Throw random cone
669  Double_t tmpRandConeEta = fJetsCont->GetJetEtaMin() + fRandom->Rndm()*TMath::Abs(fJetsCont->GetJetEtaMax()-fJetsCont->GetJetEtaMin());
670  Double_t tmpRandConePhi = fRandom->Rndm()*TMath::TwoPi();
671  Double_t tmpRandConePt = 0; // to be determined
672  Double_t tmpRandConePt3GeV = 0; // to be determined
673 
674  fAcceptedTracks = 0;
675  fTracksCont->ResetCurrentID();
676  Int_t trackcount = 0;
677  while(AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle()))
678  {
679  // Track plots
680  FillHistogramsTracks(track);
681 
682  if(IsTrackInCone(track, tmpRandConeEta, tmpRandConePhi, fJetsCont->GetJetRadius()))
683  {
684  tmpRandConePt += track->Pt();
685  if (track->Pt() > 3.0)
686  tmpRandConePt3GeV += track->Pt();
687  }
688 
689  // Add track to output array
690  trackcount++;
691  AddTrackToOutputArray(track);
692  }
693 
694  // ####### Embedding plots
695  if( (fJetOutputMode == 4) || (fJetOutputMode == 5))
696  {
697  for(Int_t i=0; i<fMatchedJets.size(); i++)
698  {
699  Double_t deltaEta = (fMatchedJets[i]->Eta()-fMatchedJetsReference[i]->Eta());
700  Double_t deltaPhi = TMath::Min(TMath::Abs(fMatchedJets[i]->Phi()-fMatchedJetsReference[i]->Phi()),TMath::TwoPi() - TMath::Abs(fMatchedJets[i]->Phi()-fMatchedJetsReference[i]->Phi()));
701  if(fMatchedJets[i]->Phi() < fMatchedJetsReference[i]->Phi())
702  deltaPhi = -deltaPhi;
703 
704  Double_t deltaR = TMath::Sqrt(deltaEta*deltaEta + deltaPhi*deltaPhi);
705 
706  FillHistogram("hEmbeddingDeltaR", fMatchedJets[i]->Pt() - fJetsCont->GetRhoVal()*fMatchedJets[i]->Area(), deltaR);
707  FillHistogram("hEmbeddingDeltaEta", fMatchedJets[i]->Pt() - fJetsCont->GetRhoVal()*fMatchedJets[i]->Area(), deltaPhi);
708  FillHistogram("hEmbeddingDeltaPhi", fMatchedJets[i]->Pt() - fJetsCont->GetRhoVal()*fMatchedJets[i]->Area(), deltaEta);
709  FillHistogram("hEmbeddingPtCorr", fMatchedJetsReference[i]->Pt(), fMatchedJets[i]->Pt() - fJetsCont->GetRhoVal()*fMatchedJets[i]->Area());
710  FillHistogram("hEmbeddingJetPt", fMatchedJetsReference[i]->Pt());
711  FillHistogram("hEmbeddingJetPhiEta", fMatchedJetsReference[i]->Phi(), fMatchedJetsReference[i]->Eta());
712  }
713  }
714 
715  // Add event to output tree
717  AddEventToTree();
718 
719  // ####### Event properties
720  FillHistogram("hRandomConePt", tmpRandConePt - fJetsCont->GetRhoVal()*fJetsCont->GetJetRadius()*fJetsCont->GetJetRadius()*TMath::Pi(), fCent);
721  FillHistogram("hRandomConePtCut3GeV", tmpRandConePt3GeV - fJetsCont->GetRhoVal()*fJetsCont->GetJetRadius()*fJetsCont->GetJetRadius()*TMath::Pi(), fCent);
722  FillHistogram("hRandomConeRawPt", tmpRandConePt, fCent);
723  FillHistogram("hRandomConeRawPtCut3GeV", tmpRandConePt3GeV, fCent);
724 
725  FillHistogram("hBackgroundPt", fJetsCont->GetRhoVal(), fCent);
726  FillHistogram("hJetCount", fAcceptedJets, fCent);
727  FillHistogram("hTrackCount", trackcount, fCent);
728  // NOTE: It is possible to use fTracksCont->GetLeadingParticle() since we do not apply additional track cuts
729  AliVTrack* leadTrack = static_cast<AliVTrack*>(fTracksCont->GetLeadingParticle());
730  if(leadTrack)
731  {
732  FillHistogram("hLeadingTrackPt", leadTrack->Pt(), fCent);
733  FillHistogram("hLeadingTrackPhi", leadTrack->Phi(), fCent);
734  FillHistogram("hLeadingTrackEta", leadTrack->Eta(), fCent);
735  FillHistogram("hLeadingTrackPhiEta", leadTrack->Phi(), leadTrack->Eta());
736  }
737 
738  return kTRUE;
739 }
740 
741 //########################################################################
742 // HELPERS
743 //########################################################################
744 
745 //________________________________________________________________________
747 {
748  if(!fPythiaInfo)
749  AliError("fPythiaInfo object not available. Is it activated with SetGeneratePythiaInfoObject()?");
750 
751  Double_t bestMatchDeltaR1 = 999.;
752  Double_t bestMatchDeltaR2 = 999.;
753 
754  fJetsCont->ResetCurrentID();
755  while(AliEmcalJet *jet = fJetsCont->GetNextAcceptJet())
756  {
757  // Check via geometrical matching if jet is connected to the initial collision
758  Double_t deltaEta1 = TMath::Abs(jet->Eta()-fPythiaInfo->GetPartonEta6());
759  Double_t deltaEta2 = TMath::Abs(jet->Eta()-fPythiaInfo->GetPartonEta7());
760  Double_t deltaPhi1 = TMath::Min(TMath::Abs(jet->Phi()-fPythiaInfo->GetPartonPhi6()),TMath::TwoPi() - TMath::Abs(jet->Phi()-fPythiaInfo->GetPartonPhi6()));
761  Double_t deltaPhi2 = TMath::Min(TMath::Abs(jet->Phi()-fPythiaInfo->GetPartonPhi7()),TMath::TwoPi() - TMath::Abs(jet->Phi()-fPythiaInfo->GetPartonPhi7()));
762 
763  Double_t deltaR1 = TMath::Sqrt(deltaEta1*deltaEta1 + deltaPhi1*deltaPhi1);
764  Double_t deltaR2 = TMath::Sqrt(deltaEta2*deltaEta2 + deltaPhi2*deltaPhi2);
765 
766  if(deltaR1 < bestMatchDeltaR1)
767  {
768  bestMatchDeltaR1 = deltaR1;
770  }
771  if(deltaR2 < bestMatchDeltaR2)
772  {
773  bestMatchDeltaR2 = deltaR2;
775  }
776  }
777 
778  if(bestMatchDeltaR1 > 0.3)
780  if(bestMatchDeltaR2 > 0.3)
782 }
783 
784 //________________________________________________________________________
786 {
787  fMatchedJets.clear();
788  fMatchedJetsReference.clear();
789 
790  // Check for a jet veto here
791  if(fJetVetoArray)
792  {
793  for(Int_t i=0; i<fJetVetoArray->GetEntries(); i++)
794  {
795  AliEmcalJet* vetoJet = static_cast<AliEmcalJet*>(fJetVetoArray->At(i));
796  UInt_t dummy = 0;
797  if(!fJetsCont->AcceptJet(vetoJet , dummy))
798  continue;
799  // if veto jet found -> return
800  if((vetoJet->Pt() >= fJetVetoMinPt) && (vetoJet->Pt() < fJetVetoMaxPt))
801  return;
802  }
803  }
804 
805  // Search for all matches above a certain threshold
806  for(Int_t i=0; i<fJetMatchingArray->GetEntries(); i++)
807  {
808  AliEmcalJet* probeJet = static_cast<AliEmcalJet*>(fJetMatchingArray->At(i));
809  UInt_t dummy = 0;
810  if(!fJetsCont->AcceptJet(probeJet , dummy))
811  continue;
812  if((probeJet->Pt() < fJetMatchingMinPt) || (probeJet->Pt() >= fJetMatchingMaxPt))
813  continue;
814 
815  AliEmcalJet* matchedJet = 0;
816  AliEmcalJet* matchedJetReference = 0;
817  Double_t bestMatchDeltaR = 999.;
818  fJetsCont->ResetCurrentID();
819  // Loop over all embedded jets to find the best match
820  while(AliEmcalJet* embeddedJet = fJetsCont->GetNextAcceptJet())
821  {
822  Double_t deltaEta = (embeddedJet->Eta()-probeJet->Eta());
823  Double_t deltaPhi = TMath::Min(TMath::Abs(embeddedJet->Phi()-probeJet->Phi()),TMath::TwoPi() - TMath::Abs(embeddedJet->Phi()-probeJet->Phi()));
824  Double_t deltaR = TMath::Sqrt(deltaEta*deltaEta + deltaPhi*deltaPhi);
825 
826  // Cut jets too far away
827  if (deltaR > fJetMatchingMaxDistance)
828  continue;
829  // Cut jets with too small pT
830  if(embeddedJet->Pt() - fJetsCont->GetRhoVal()*embeddedJet->Area() < fJetMatchingMinSharedFraction*probeJet->Pt())
831  continue;
832  // Cut jets with too high pT
833  if(embeddedJet->Pt() - fJetsCont->GetRhoVal()*embeddedJet->Area() > (fJetMatchingMaxSharedFraction*probeJet->Pt() + fJetMatchingMaxEmbeddingOffset))
834  continue;
835 
836  // Search for the best match
837  if(deltaR < bestMatchDeltaR)
838  {
839  bestMatchDeltaR = deltaR;
840  matchedJet = embeddedJet;
841  matchedJetReference = probeJet;
842  }
843  }
844  // Put matched jet to a list
845  if(matchedJet && matchedJetReference)
846  {
847  fMatchedJets.push_back(matchedJet);
848  fMatchedJetsReference.push_back(matchedJetReference);
849  }
850  }
851 
852  // ############ On demand, search for matches of leading and subleading jet in MC and delete rest
854  {
855  Int_t jetLeadingIndex = -1;
856  Int_t jetSubLeadingIndex = -1;
857  Double_t tmpLeadingPt = 0;
858  Double_t tmpSubleadingPt = 0;
859 
860  // Search leading/subleading in matched MC jets
861  for(Int_t i=0; i<fMatchedJetsReference.size(); i++)
862  {
863  AliEmcalJet* matchedJet = fMatchedJetsReference[i];
864  if (matchedJet->Pt() > tmpLeadingPt)
865  {
866  jetSubLeadingIndex = jetLeadingIndex;
867  jetLeadingIndex = i;
868  tmpSubleadingPt = tmpLeadingPt;
869  tmpLeadingPt = matchedJet->Pt();
870  }
871  else if (matchedJet->Pt() > tmpSubleadingPt)
872  {
873  jetSubLeadingIndex = i;
874  tmpSubleadingPt = matchedJet->Pt();
875  }
876  }
877 
878  AliEmcalJet* matchedLeading = 0;
879  AliEmcalJet* refLeading = 0;
880  AliEmcalJet* matchedSubLeading = 0;
881  AliEmcalJet* refSubLeading = 0;
882 
883  // Cache leading jet, if found
884  if(jetLeadingIndex >= 0)
885  {
886  matchedLeading = fMatchedJets[jetLeadingIndex];
887  refLeading = fMatchedJetsReference[jetLeadingIndex];
888  }
889  // Cache subleading jet, if found
890  if(jetSubLeadingIndex >= 0)
891  {
892  matchedSubLeading = fMatchedJets[jetSubLeadingIndex];
893  refSubLeading = fMatchedJetsReference[jetSubLeadingIndex];
894  }
895 
896  // Delete old list
897  fMatchedJets.clear();
898  fMatchedJetsReference.clear();
899 
900  // Add jets
901  if(matchedLeading)
902  {
903  fMatchedJets.push_back(matchedLeading);
904  fMatchedJetsReference.push_back(refLeading);
905  }
906  if(matchedSubLeading && fJetMatchingUseOnlyNLeading >= 2)
907  {
908  fMatchedJets.push_back(matchedSubLeading);
909  fMatchedJetsReference.push_back(refSubLeading);
910  }
911  }
912 
913 }
914 
915 //________________________________________________________________________
917 {
918  // This is to use a full cone in phi even at the edges of phi (2pi -> 0) (0 -> 2pi)
919  Double_t trackPhi = 0.0;
920  if (track->Phi() > (TMath::TwoPi() - (radius-phi)))
921  trackPhi = track->Phi() - TMath::TwoPi();
922  else if (track->Phi() < (phi+radius - TMath::TwoPi()))
923  trackPhi = track->Phi() + TMath::TwoPi();
924  else
925  trackPhi = track->Phi();
926 
927  if ( TMath::Abs(trackPhi-phi)*TMath::Abs(trackPhi-phi) + TMath::Abs(track->Eta()-eta)*TMath::Abs(track->Eta()-eta) <= radius*radius)
928  return kTRUE;
929 
930  return kFALSE;
931 }
932 
933 
934 //________________________________________________________________________
936 {
937  // Calculate leading + subleading jet
939  if(fJetOutputMode==6)
941  else if(fJetOutputMode==4 || fJetOutputMode==5)
942  GetMatchingJets();
943 }
944 
945 //________________________________________________________________________
947 {
948  Double_t fakeFactor = 0;
949 
950  // Loop over all jet constituents
951  for(Int_t i = 0; i < jet->GetNumberOfTracks(); i++)
952  {
953  AliVParticle* constituent = static_cast<AliVParticle*>(jet->TrackAt(i, fTracksCont->GetArray()));
954 
955  Double_t deltaPhi = TMath::Min(TMath::Abs(jet->Phi()-constituent->Phi()),TMath::TwoPi() - TMath::Abs(jet->Phi()-constituent->Phi()));
956  Double_t deltaR = TMath::Sqrt( (jet->Eta() - constituent->Eta())*(jet->Eta() - constituent->Eta()) + deltaPhi*deltaPhi );
957  fakeFactor += constituent->Pt() * TMath::Sin(deltaR);
958  }
959 
960  return fakeFactor;
961 }
962 
963 //________________________________________________________________________
965 {
966  fEventCriteriumMode = type;
967 
968  if(fEventCriteriumMode==0)
969  AliWarning("Set event criterium to 'default' -- no further selection criterium.");
970  else if(fEventCriteriumMode==1)
971  AliWarning("Set event criterium to 'background' -- select events with certain backgrounds");
972  else if(fEventCriteriumMode==2)
973  AliWarning("Set event criterium to 'simple jet trigger' -- select events with certain minimum leading jet pT (bgrd corr.)");
974  else if(fEventCriteriumMode==3)
975  AliWarning("Set event criterium to 'simple dijet trigger' -- select events with certain minimum leading + subleading jet pT (bgrd corr.)");
976  else
977  {
978  AliFatal("Event criterium not valid.");
979  }
980 }
981 
982 
983 //________________________________________________________________________
984 void AliAnalysisTaskChargedJetsHadronCF::GetLeadingJets(const char* opt, AliEmcalJet*& jetLeading, AliEmcalJet*& jetSubLeading)
985 {
986  // Customized from AliJetContainer::GetLeadingJet()
987  // Get the leading+subleading jet; if opt contains "rho" the sorting is according to pt-A*rho
988 
989  TString option(opt);
990  option.ToLower();
991 
992  jetLeading = 0;
993  jetSubLeading = 0;
994 
995  fJetsCont->ResetCurrentID();
996  Double_t tmpLeadingPt = 0;
997  Double_t tmpSubleadingPt = 0;
998 
999  if (option.Contains("rho")) {
1000  while (AliEmcalJet* jet = fJetsCont->GetNextAcceptJet()) {
1001  if ( (jet->Pt()-jet->Area()*fJetsCont->GetRhoVal()) > tmpLeadingPt )
1002  {
1003  jetSubLeading = jetLeading;
1004  jetLeading = jet;
1005  tmpSubleadingPt = tmpLeadingPt;
1006  tmpLeadingPt = jet->Pt()-jet->Area()*fJetsCont->GetRhoVal();
1007  }
1008  else if ( (jet->Pt()-jet->Area()*fJetsCont->GetRhoVal()) > tmpSubleadingPt )
1009  {
1010  jetSubLeading = jet;
1011  tmpSubleadingPt = jet->Pt()-jet->Area()*fJetsCont->GetRhoVal();
1012  }
1013  }
1014  }
1015  else {
1016  while (AliEmcalJet* jet = fJetsCont->GetNextAcceptJet()) {
1017  if ( (jet->Pt()) > tmpLeadingPt )
1018  {
1019  jetSubLeading = jetLeading;
1020  jetLeading = jet;
1021  tmpSubleadingPt = tmpLeadingPt;
1022  tmpLeadingPt = jet->Pt();
1023  }
1024  else if ( (jet->Pt()) > tmpSubleadingPt )
1025  {
1026  jetSubLeading = jet;
1027  tmpSubleadingPt = jet->Pt();
1028  }
1029  }
1030  }
1031 }
1032 
1033 //________________________________________________________________________
1035 {
1036  // Method for the correct logarithmic binning of histograms
1037  TAxis *axis = h->GetAxis(axisNumber);
1038  int bins = axis->GetNbins();
1039 
1040  Double_t from = axis->GetXmin();
1041  Double_t to = axis->GetXmax();
1042  Double_t *newBins = new Double_t[bins + 1];
1043 
1044  newBins[0] = from;
1045  Double_t factor = pow(to/from, 1./bins);
1046 
1047  for (int i = 1; i <= bins; i++) {
1048  newBins[i] = factor * newBins[i-1];
1049  }
1050  axis->Set(bins, newBins);
1051  delete [] newBins;
1052 }
1053 
1054 //________________________________________________________________________
1056 {
1057  TH1* tmpHist = static_cast<TH1*>(fOutput->FindObject(key));
1058  if(!tmpHist)
1059  {
1060  AliError(Form("Cannot find histogram <%s> ",key)) ;
1061  return;
1062  }
1063 
1064  tmpHist->Fill(x);
1065 }
1066 
1067 //________________________________________________________________________
1069 {
1070  TH1* tmpHist = static_cast<TH1*>(fOutput->FindObject(key));
1071  if(!tmpHist)
1072  {
1073  AliError(Form("Cannot find histogram <%s> ",key));
1074  return;
1075  }
1076 
1077  if (tmpHist->IsA()->GetBaseClass("TH1"))
1078  static_cast<TH1*>(tmpHist)->Fill(x,y); // Fill x with y
1079  else if (tmpHist->IsA()->GetBaseClass("TH2"))
1080  static_cast<TH2*>(tmpHist)->Fill(x,y); // Fill x,y with 1
1081 }
1082 
1083 //________________________________________________________________________
1085 {
1086  TH2* tmpHist = static_cast<TH2*>(fOutput->FindObject(key));
1087  if(!tmpHist)
1088  {
1089  AliError(Form("Cannot find histogram <%s> ",key));
1090  return;
1091  }
1092 
1093  tmpHist->Fill(x,y,add);
1094 }
1095 
1096 //________________________________________________________________________
1098 {
1099  TH3* tmpHist = static_cast<TH3*>(fOutput->FindObject(key));
1100  if(!tmpHist)
1101  {
1102  AliError(Form("Cannot find histogram <%s> ",key));
1103  return;
1104  }
1105 
1106  if(add)
1107  tmpHist->Fill(x,y,z,add);
1108  else
1109  tmpHist->Fill(x,y,z);
1110 }
1111 
1112 
1113 //________________________________________________________________________
1114 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)
1115 {
1116  T* tmpHist = new T(name, title, xBins, xMin, xMax);
1117 
1118  tmpHist->GetXaxis()->SetTitle(xTitle);
1119  tmpHist->GetYaxis()->SetTitle(yTitle);
1120  tmpHist->SetOption(options);
1121  tmpHist->SetMarkerStyle(kFullCircle);
1122  tmpHist->Sumw2();
1123 
1124  fOutput->Add(tmpHist);
1125 
1126  return tmpHist;
1127 }
1128 
1129 //________________________________________________________________________
1130 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)
1131 {
1132  T* tmpHist = new T(name, title, xBins, xMin, xMax, yBins, yMin, yMax);
1133  tmpHist->GetXaxis()->SetTitle(xTitle);
1134  tmpHist->GetYaxis()->SetTitle(yTitle);
1135  tmpHist->GetZaxis()->SetTitle(zTitle);
1136  tmpHist->SetOption(options);
1137  tmpHist->SetMarkerStyle(kFullCircle);
1138  tmpHist->Sumw2();
1139 
1140  fOutput->Add(tmpHist);
1141 
1142  return tmpHist;
1143 }
1144 
1145 //________________________________________________________________________
1146 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)
1147 {
1148  T* tmpHist = new T(name, title, xBins, xMin, xMax, yBins, yMin, yMax, zBins, zMin, zMax);
1149  tmpHist->GetXaxis()->SetTitle(xTitle);
1150  tmpHist->GetYaxis()->SetTitle(yTitle);
1151  tmpHist->GetZaxis()->SetTitle(zTitle);
1152  tmpHist->SetOption(options);
1153  tmpHist->SetMarkerStyle(kFullCircle);
1154  tmpHist->Sumw2();
1155 
1156  fOutput->Add(tmpHist);
1157 
1158  return tmpHist;
1159 }
1160 
1161 //________________________________________________________________________
1163 {
1164  // Called once at the end of the analysis.
1165 }
1166 
Double_t fEventExtractionMaxJetPt
maximum jet pt of recorded events
Short_t Charge() const
Definition: AliEmcalJet.h:110
Double_t fEventCriteriumMinSubleadingJetPt
Min subleading jet.
Double_t Area() const
Definition: AliEmcalJet.h:117
Double_t GetRhoVal() const
double Double_t
Definition: External.C:58
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
TString fJetMatchingArrayName
Name of array used to match jets.
AliJetContainer * GetJetContainer(Int_t i=0) const
Definition: External.C:244
Float_t GetPartonEta7() const
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
Double_t fJetMatchingMaxEmbeddingOffset
An embedded jet must NOT carry more than max fraction + max embedding offset.
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
Simple class containing basic information for a constituent.
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
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:126
UShort_t T(UShort_t m, UShort_t t)
Definition: RingBits.C:60
Double_t fEventCriteriumMinJetDeltaPhi
Min jet delta phi in dijet criterium.
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Int_t GetPartonFlag6() const
void AddJetConstituent(Float_t eta, Float_t phi, Float_t pt, Short_t charge, Short_t pid=0)
AliParticleContainer * GetParticleContainer() const
int Int_t
Definition: External.C:63
Double_t fJetMatchingMaxPt
Max pt cut applied on the matchArray jets.
unsigned int UInt_t
Definition: External.C:33
TString fJetParticleArrayName
Name of fJetsOutput array.
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")
Int_t fNumberOfCentralityBins
Number of centrality bins.
virtual AliVParticle * GetLeadingParticle(const char *opt="")
std::vector< AliEmcalJet * > fMatchedJets
Jets matched in an event (embedded)
Double_t fExtractionPercentage
percentage that is recorded
Double_t fCent
!event centrality
Double_t fJetVetoMinPt
Min pt cut applied on the fJetVetoArray jets.
Double_t fJetVetoMaxPt
Max pt cut applied on the fJetVetoArray jets.
TClonesArray * fJetsOutput
! Array of basic correlation particles attached to the event (jets)
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()
TTree * fJetsTree
! Jets that will be saved to a tree (optionally)
Double_t fJetMatchingMaxSharedFraction
An embedded jet must NOT carry more than max fraction + max embedding offset.
Double_t Pt() const
Definition: AliEmcalJet.h:96
Bool_t IsTrackInCone(AliVParticle *track, Double_t eta, Double_t phi, Double_t radius)
TString fJetVetoArrayName
Name of array used for veto jets.
Float_t GetJetRadius() const
AliEmcalList * fOutput
!output list
TClonesArray * fJetMatchingArray
! Array of generated jets imported into task (toy model/embedding)
Short_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:147
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.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:44
Declaration of class AliEmcalPythiaInfo.
const char Option_t
Definition: External.C:48
Double_t fExtractionMaxPt
maximum pt of recorded jets
Double_t fJetMatchingMaxDistance
Max distance allowed to accept a matching jet (for embedding)
virtual AliVParticle * GetNextAcceptParticle()
bool Bool_t
Definition: External.C:53
Double_t fJetMatchingMinPt
Min pt cut applied on the matchArray jets.
Int_t fJetMatchingUseOnlyNLeading
Number of matched leading jets that will be used.
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)
Double_t fJetMatchingMinSharedFraction
An embedded jet must carry this pt fraction to be accepted.