AliPhysics  master (3d17d9d)
AliAnalysisTaskEmcalSoftDropResponse.cxx
Go to the documentation of this file.
1 /************************************************************************************
2  * Copyright (C) 2019, Copyright Holders of the ALICE Collaboration *
3  * All rights reserved. *
4  * *
5  * Redistribution and use in source and binary forms, with or without *
6  * modification, are permitted provided that the following conditions are met: *
7  * * Redistributions of source code must retain the above copyright *
8  * notice, this list of conditions and the following disclaimer. *
9  * * Redistributions in binary form must reproduce the above copyright *
10  * notice, this list of conditions and the following disclaimer in the *
11  * documentation and/or other materials provided with the distribution. *
12  * * Neither the name of the <organization> nor the *
13  * names of its contributors may be used to endorse or promote products *
14  * derived from this software without specific prior written permission. *
15  * *
16  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND *
17  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED *
18  * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE *
19  * DISCLAIMED. IN NO EVENT SHALL ALICE COLLABORATION BE LIABLE FOR ANY *
20  * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES *
21  * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; *
22  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND *
23  * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT *
24  * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS *
25  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
26  ************************************************************************************/
27 #include <iostream>
28 #include <memory>
29 #include <vector>
30 
31 #include <TArrayD.h>
32 #include <TBinning.h>
33 #include <TCustomBinning.h>
34 #include <TH1.h>
35 #include <TH2.h>
36 #include <TLinearBinning.h>
37 #include <TLorentzVector.h>
38 #include <TRandom.h>
39 
40 #include <fastjet/ClusterSequence.hh>
41 #include <fastjet/PseudoJet.hh>
42 #include <fastjet/contrib/SoftDrop.hh>
43 #include <fastjet/config.h>
44 #if FASJET_VERSION_NUMBER >= 30302
45 #include <fastjet/tools/Recluster.hh>
46 #else
47 #include <fastjet/contrib/Recluster.hh>
48 #endif
49 
50 #include <RooUnfoldResponse.h>
51 
52 #include <AliAODEvent.h>
53 #include <AliAODInputHandler.h>
54 #include <AliAnalysisManager.h>
56 #include <AliClusterContainer.h>
57 #include <AliEmcalJet.h>
59 #include <AliJetContainer.h>
60 #include <AliMCParticleContainer.h>
61 #include <AliTrackContainer.h>
62 #include <AliVCluster.h>
63 #include <AliVParticle.h>
64 #include <AliVTrack.h>
65 
67 
68  using namespace PWGJE::EMCALJetTasks;
69 
71  fBinningMode(kSDModeINT7),
72  fFractionResponseClosure(0.5),
73  fZcut(0.1),
74  fBeta(0.),
75  fReclusterizer(kCAAlgo),
76  fSampleFraction(1.),
77  fMinFractionShared(0),
78  fHasResponseMatrixSparse(false),
79  fHasResponseMatrixRooUnfold(true),
80  fUseChargedConstituents(true),
81  fUseNeutralConstituents(true),
82  fNameMCParticles("mcparticles"),
83  fSampleSplitter(nullptr),
84  fSampleTrimmer(nullptr),
85  fPartLevelPtBinning(nullptr),
86  fDetLevelPtBinning(nullptr),
87  fIsEmbeddedEvent(false),
88  fNamePartLevelJetContainer(""),
89  fNameDetLevelJetContainer(""),
90  fNameUnSubLevelJetContainer(""),
91  fZgResponse(),
92  fZgResponseClosure(),
93  fRgResponse(),
94  fRgResponseClosure(),
95  fNsdResponse(),
96  fNsdResponseClosure(),
97  fThetagResponse(),
98  fThetagResponseClosure(),
99  fHistManager("AliAnalysisTaskSoftDropResponse")
100 {
101 }
102 
103 AliAnalysisTaskEmcalSoftDropResponse::AliAnalysisTaskEmcalSoftDropResponse(const char *name) : AliAnalysisTaskEmcalJet(name, kTRUE),
104  fBinningMode(kSDModeINT7),
105  fFractionResponseClosure(0.5),
106  fZcut(0.1),
107  fBeta(0.),
108  fReclusterizer(kCAAlgo),
109  fSampleFraction(1.),
110  fMinFractionShared(0),
111  fHasResponseMatrixSparse(false),
112  fHasResponseMatrixRooUnfold(true),
113  fUseChargedConstituents(true),
114  fUseNeutralConstituents(true),
115  fNameMCParticles("mcparticles"),
116  fSampleSplitter(nullptr),
117  fSampleTrimmer(nullptr),
118  fPartLevelPtBinning(nullptr),
119  fDetLevelPtBinning(nullptr),
120  fNamePartLevelJetContainer(""),
121  fNameDetLevelJetContainer(""),
122  fNameUnSubLevelJetContainer(""),
123  fIsEmbeddedEvent(false),
124  fZgResponse(),
125  fZgResponseClosure(),
126  fRgResponse(),
127  fRgResponseClosure(),
128  fNsdResponse(),
129  fNsdResponseClosure(),
130  fThetagResponse(),
131  fThetagResponseClosure(),
132  fHistManager(name)
133 {
135 }
136 
138 {
140  delete fPartLevelPtBinning;
141  if (fDetLevelPtBinning)
142  delete fDetLevelPtBinning;
143  if (fSampleSplitter)
144  delete fSampleSplitter;
145  if (fSampleTrimmer)
146  delete fSampleTrimmer;
147 }
148 
150 {
153 
154  fSampleSplitter = new TRandom;
155  if (fSampleFraction < 1.)
156  fSampleTrimmer = new TRandom;
157 
158  if (!fPartLevelPtBinning)
160  if (!fDetLevelPtBinning)
162  std::unique_ptr<TBinning> zgbinning(GetZgBinning()),
163  rgbinning(GetRgBinning(R)),
164  nsdbinning(new TLinearBinning(22, -1.5, 20.5)), // Negative bins are for untagged jets
165  thetagbinning(new TLinearBinning(11, -0.1, 1.)),
166  ptbinningFine(new TLinearBinning(500, 0., 500.));
167  TArrayD binEdgesZg, binEdgesRg, binEdgesNsd, binEdgesThetag, binEdgesPtPart, binEdgesPtDet, binEdgesPtFine;
168  zgbinning->CreateBinEdges(binEdgesZg);
169  rgbinning->CreateBinEdges(binEdgesRg);
170  nsdbinning->CreateBinEdges(binEdgesNsd);
171  thetagbinning->CreateBinEdges(binEdgesThetag);
172  ptbinningFine->CreateBinEdges(binEdgesPtFine);
173  fPartLevelPtBinning->CreateBinEdges(binEdgesPtPart);
174  fDetLevelPtBinning->CreateBinEdges(binEdgesPtDet);
175 
176  const TBinning *sparsebinningZg[4] = {zgbinning.get(), ptbinningFine.get(), zgbinning.get(), ptbinningFine.get()},
177  *sparsebinningRg[4] = {rgbinning.get(), ptbinningFine.get(), rgbinning.get(), ptbinningFine.get()},
178  *sparsebinningNsd[4] = {nsdbinning.get(), ptbinningFine.get(), nsdbinning.get(), ptbinningFine.get()},
179  *sparsebinningThetag[4] = {thetagbinning.get(), ptbinningFine.get(), thetagbinning.get(), ptbinningFine.get()};
180 
181  //Need to do centrality bins in the histograms if it is not pp
182  if (fForceBeamType != kpp)
183  {
184  for (Int_t cent = 0; cent < fNcentBins; cent++)
185  {
187  fHistManager.CreateTH2(Form("hZgDetLevel_%d", cent), Form("Zg response at detector level, %d centrality bin", cent), binEdgesZg.GetSize() - 1, binEdgesZg.GetArray(), binEdgesPtDet.GetSize() - 1, binEdgesPtDet.GetArray());
188  fHistManager.CreateTH2(Form("hZgPartLevel_%d", cent), Form("Zg response at particle level, %d centrality bin", cent), binEdgesZg.GetSize() - 1, binEdgesZg.GetArray(), binEdgesPtPart.GetSize() - 1, binEdgesPtPart.GetArray());
189  fHistManager.CreateTH2(Form("hZgPartLevelTruncated_%d", cent), Form("Zg response at particle level (truncated), %d centrality bin", cent), binEdgesZg.GetSize() - 1, binEdgesZg.GetArray(), binEdgesPtPart.GetSize() - 1, binEdgesPtPart.GetArray());
190  fHistManager.CreateTH2(Form("hRgDetLevel_%d", cent), Form("Rg response at detector level, %d centrality bin", cent), binEdgesRg.GetSize() - 1, binEdgesRg.GetArray(), binEdgesPtDet.GetSize() - 1, binEdgesPtDet.GetArray());
191  fHistManager.CreateTH2(Form("hRgPartLevel_%d", cent), Form("Rg response at particle level, %d centrality bin", cent), binEdgesRg.GetSize() - 1, binEdgesRg.GetArray(), binEdgesPtPart.GetSize() - 1, binEdgesPtPart.GetArray());
192  fHistManager.CreateTH2(Form("hRgPartLevelTruncated_%d", cent), Form("Rg response at particle level (truncated), %d centrality bin", cent), binEdgesRg.GetSize() - 1, binEdgesRg.GetArray(), binEdgesPtPart.GetSize() - 1, binEdgesPtPart.GetArray());
193  fHistManager.CreateTH2(Form("hNsdDetLevel_%d", cent), Form("Nsd response at detector level, %d centrality bin", cent), binEdgesNsd.GetSize() - 1, binEdgesNsd.GetArray(), binEdgesPtDet.GetSize() - 1, binEdgesPtDet.GetArray());
194  fHistManager.CreateTH2(Form("hNsdPartLevel_%d", cent), Form("Nsd response at particle level, %d centrality bin", cent), binEdgesNsd.GetSize() - 1, binEdgesNsd.GetArray(), binEdgesPtPart.GetSize() - 1, binEdgesPtPart.GetArray());
195  fHistManager.CreateTH2(Form("hNsdPartLevelTruncated_%d", cent), Form("Nsd response at particle level (truncated), %d centrality bin", cent), binEdgesNsd.GetSize() - 1, binEdgesNsd.GetArray(), binEdgesPtPart.GetSize() - 1, binEdgesPtPart.GetArray());
196  fHistManager.CreateTH2(Form("hThetagDetLevel_%d", cent), Form("Thetag response at detector level, %d centrality bin", cent), binEdgesThetag.GetSize() - 1, binEdgesThetag.GetArray(), binEdgesPtDet.GetSize() - 1, binEdgesPtDet.GetArray());
197  fHistManager.CreateTH2(Form("hThetagPartLevel_%d", cent), Form("Thetag response at particle level, %d centrality bin", cent), binEdgesThetag.GetSize() - 1, binEdgesThetag.GetArray(), binEdgesPtPart.GetSize() - 1, binEdgesPtPart.GetArray());
198  fHistManager.CreateTH2(Form("hThetagPartLevelTruncated_%d", cent), Form("Thetag response at particle level (truncated), %d centrality bin", cent), binEdgesThetag.GetSize() - 1, binEdgesThetag.GetArray(), binEdgesPtPart.GetSize() - 1, binEdgesPtPart.GetArray());
199 
200  // For closure test
201  fHistManager.CreateTH2(Form("hZgPartLevelClosureNoResp_%d", cent), Form("Zg response at particle level (closure test, jets not used for the response matrix), %d centrality bin", cent), binEdgesZg.GetSize() - 1, binEdgesZg.GetArray(), binEdgesPtPart.GetSize() - 1, binEdgesPtPart.GetArray());
202  fHistManager.CreateTH2(Form("hZgDetLevelClosureNoResp_%d", cent), Form("Zg response at detector level (closure test, jets not used for the response matrix), %d centrality bin", cent), binEdgesZg.GetSize() - 1, binEdgesZg.GetArray(), binEdgesPtDet.GetSize() - 1, binEdgesPtDet.GetArray());
203  fHistManager.CreateTH2(Form("hZgPartLevelClosureResp_%d", cent), Form("Zg response at particle level (closure test, jets used for the response matrix), %d centrality bin", cent), binEdgesZg.GetSize() - 1, binEdgesZg.GetArray(), binEdgesPtPart.GetSize() - 1, binEdgesPtPart.GetArray());
204  fHistManager.CreateTH2(Form("hZgDetLevelClosureResp_%d", cent), Form("Zg response at detector level (closure test, jets used for the response matrix), %d centrality bin", cent), binEdgesZg.GetSize() - 1, binEdgesZg.GetArray(), binEdgesPtDet.GetSize() - 1, binEdgesPtDet.GetArray());
205  fHistManager.CreateTH2(Form("hRgPartLevelClosureNoResp_%d", cent), Form("Rg response at particle level (closure test, jets not used for the response matrix), %d centrality bin", cent), binEdgesRg.GetSize() - 1, binEdgesRg.GetArray(), binEdgesPtPart.GetSize() - 1, binEdgesPtPart.GetArray());
206  fHistManager.CreateTH2(Form("hRgDetLevelClosureNoResp_%d", cent), Form("Rg response at detector level (closure test, jets not used for the response matrix), %d centrality bin", cent), binEdgesRg.GetSize() - 1, binEdgesRg.GetArray(), binEdgesPtDet.GetSize() - 1, binEdgesPtDet.GetArray());
207  fHistManager.CreateTH2(Form("hRgPartLevelClosureResp_%d", cent), Form("Rg response at particle level (closure test, jets used for the response matrix), %d centrality bin", cent), binEdgesRg.GetSize() - 1, binEdgesRg.GetArray(), binEdgesPtPart.GetSize() - 1, binEdgesPtPart.GetArray());
208  fHistManager.CreateTH2(Form("hRgDetLevelClosureResp_%d", cent), Form("Rg response at detector level (closure test, jets used for the response matrix), %d centrality bin", cent), binEdgesRg.GetSize() - 1, binEdgesRg.GetArray(), binEdgesPtDet.GetSize() - 1, binEdgesPtDet.GetArray());
209  fHistManager.CreateTH2(Form("hNsdPartLevelClosureNoResp_%d", cent), Form("Nsd response at particle level (closure test, jets not used for the response matrix), %d centrality bin", cent), binEdgesNsd.GetSize() - 1, binEdgesNsd.GetArray(), binEdgesPtPart.GetSize() - 1, binEdgesPtPart.GetArray());
210  fHistManager.CreateTH2(Form("hNsdDetLevelClosureNoResp_%d", cent), Form("Nsd response at detector level (closure test, jets not used for the response matrix), %d centrality bin", cent), binEdgesNsd.GetSize() - 1, binEdgesNsd.GetArray(), binEdgesPtDet.GetSize() - 1, binEdgesPtDet.GetArray());
211  fHistManager.CreateTH2(Form("hNsdPartLevelClosureResp_%d", cent), Form("Nsd response at particle level (closure test, jets used for the response matrix), %d centrality bin", cent), binEdgesNsd.GetSize() - 1, binEdgesNsd.GetArray(), binEdgesPtPart.GetSize() - 1, binEdgesPtPart.GetArray());
212  fHistManager.CreateTH2(Form("hNsdDetLevelClosureResp_%d", cent), Form("Nsd response at detector level (closure test, jets used for the response matrix), %d centrality bin", cent), binEdgesNsd.GetSize() - 1, binEdgesNsd.GetArray(), binEdgesPtDet.GetSize() - 1, binEdgesPtDet.GetArray());
213  fHistManager.CreateTH2(Form("hThetagPartLevelClosureNoResp_%d", cent), Form("Thetag response at particle level (closure test, jets not used for the response matrix), %d centrality bin", cent), binEdgesThetag.GetSize() - 1, binEdgesThetag.GetArray(), binEdgesPtPart.GetSize() - 1, binEdgesPtPart.GetArray());
214  fHistManager.CreateTH2(Form("hThetagDetLevelClosureNoResp_%d", cent), Form("Thetag response at detector level (closure test, jets not used for the response matrix), %d centrality bin", cent), binEdgesThetag.GetSize() - 1, binEdgesThetag.GetArray(), binEdgesPtDet.GetSize() - 1, binEdgesPtDet.GetArray());
215  fHistManager.CreateTH2(Form("hThetagPartLevelClosureResp_%d", cent), Form("Thetag response at particle level (closure test, jets used for the response matrix), %d centrality bin", cent), binEdgesThetag.GetSize() - 1, binEdgesThetag.GetArray(), binEdgesPtPart.GetSize() - 1, binEdgesPtPart.GetArray());
216  fHistManager.CreateTH2(Form("hThetagDetLevelClosureResp_%d", cent), Form("Thetag response at detector level (closure test, jets used for the response matrix), %d centrality bin", cent), binEdgesThetag.GetSize() - 1, binEdgesThetag.GetArray(), binEdgesPtDet.GetSize() - 1, binEdgesPtDet.GetArray());
217 
218  RooUnfoldResponse *r_zg = new RooUnfoldResponse(Form("hZgResponse_%d", cent), Form("z_{g} response matrix, %d centrality bin", cent)),
219  *r_rg = new RooUnfoldResponse(Form("hRgResponse_%d", cent), Form("r_{g} response matrix, %d centrality bin", cent)),
220  *r_nsd = new RooUnfoldResponse(Form("hNsdResponse_%d", cent), Form("n_{SD} response matrix, %d centrality bin", cent)),
221  *r_thetag = new RooUnfoldResponse(Form("hThetagResponse_%d", cent), Form("#Theta_{g} response matrix, %d centrality bin", cent)),
222  *r_zg_closure = new RooUnfoldResponse(Form("hZgResponseClosure_%d", cent), Form("z_{g} response matrix for the closure test, %d centrality bin", cent)),
223  *r_rg_closure = new RooUnfoldResponse(Form("hRgResponseClosure_%d", cent), Form("r_{g} response matrix for the closure test, %d centrality bin", cent)),
224  *r_nsd_closure = new RooUnfoldResponse(Form("hNsdResponseClosure_%d", cent), Form("n_{SD} response matrix for the closure test, %d centrality bin", cent)),
225  *r_thetag_closure = new RooUnfoldResponse(Form("hThetagResponseClosure_%d", cent), Form("#Theta_{g} response matrix for the closure test, %d centrality bin", cent));
226  TString nameZgDetLevel = TString::Format("hZgDetLevel_%d", cent);
227  TString nameZgPartLevel = TString::Format("hZgPartLevel_%d", cent);
228  r_zg->Setup((TH1 *)fHistManager.FindObject(nameZgDetLevel), (TH1 *)fHistManager.FindObject(nameZgPartLevel));
229  r_zg_closure->Setup((TH1 *)fHistManager.FindObject(nameZgDetLevel), (TH1 *)fHistManager.FindObject(nameZgPartLevel));
230  TString nameRgDetLevel = TString::Format("hRgDetLevel_%d", cent);
231  TString nameRgPartLevel = TString::Format("hRgPartLevel_%d", cent);
232  r_rg->Setup((TH1 *)fHistManager.FindObject(nameRgDetLevel), (TH1 *)fHistManager.FindObject(nameRgPartLevel));
233  r_rg_closure->Setup((TH1 *)fHistManager.FindObject(nameRgDetLevel), (TH1 *)fHistManager.FindObject(nameRgPartLevel));
234  TString nameNsdDetLevel = TString::Format("hNsdDetLevel_%d", cent);
235  TString nameNsdPartLevel = TString::Format("hNsdPartLevel_%d", cent);
236  r_nsd->Setup((TH1 *)fHistManager.FindObject(nameNsdDetLevel), (TH1 *)fHistManager.FindObject(nameNsdPartLevel));
237  r_nsd_closure->Setup((TH1 *)fHistManager.FindObject(nameNsdDetLevel), (TH1 *)fHistManager.FindObject(nameNsdPartLevel));
238  TString nameThetagDetLevel = TString::Format("hThetagDetLevel_%d", cent);
239  TString nameThetagPartLevel = TString::Format("hThetagPartLevel_%d", cent);
240  r_thetag->Setup((TH1 *)fHistManager.FindObject(nameThetagDetLevel), (TH1 *)fHistManager.FindObject(nameThetagPartLevel));
241  r_thetag_closure->Setup((TH1 *)fHistManager.FindObject(nameThetagDetLevel), (TH1 *)fHistManager.FindObject(nameThetagPartLevel));
242  fZgResponse.push_back(r_zg);
243  fZgResponseClosure.push_back(r_zg_closure);
244  fRgResponse.push_back(r_rg);
245  fRgResponseClosure.push_back(r_rg_closure);
246  fNsdResponse.push_back(r_nsd);
247  fNsdResponseClosure.push_back(r_nsd_closure);
248  fThetagResponse.push_back(r_thetag);
249  fThetagResponseClosure.push_back(r_thetag_closure);
250  }
252  fHistManager.CreateTHnSparse(Form("hZgResponseSparse_%d", cent), Form("z_{g} response matrix, %d centrality bin", cent), 4, sparsebinningZg);
253  fHistManager.CreateTHnSparse(Form("hRgResponseSparse_%d", cent), Form("z_{g} response matrix, %d centrality bin", cent), 4, sparsebinningRg);
254  fHistManager.CreateTHnSparse(Form("hNsdResponseSparse_%d", cent), Form("z_{g} response matrix, %d centrality bin", cent), 4, sparsebinningNsd);
255  fHistManager.CreateTHnSparse(Form("hThetagResponseSparse_%d", cent), Form("z_{g} response matrix, %d centrality bin", cent), 4, sparsebinningThetag);
256  fHistManager.CreateTHnSparse(Form("hZgResponseClosureSparse_%d", cent), Form("z_{g} response matrix for closure test, %d centrality bin", cent), 4, sparsebinningZg);
257  fHistManager.CreateTHnSparse(Form("hRgResponseClosureSparse_%d", cent), Form("z_{g} response matrix for closure test, %d centrality bin", cent), 4, sparsebinningRg);
258  fHistManager.CreateTHnSparse(Form("hNsdResponseClosureSparse_%d", cent), Form("z_{g} response matrix for closure test, %d centrality bin", cent), 4, sparsebinningNsd);
259  fHistManager.CreateTHnSparse(Form("hThetagResponseClosureSparse_%d", cent), Form("z_{g} response matrix for closure test, %d centrality bin", cent), 4, sparsebinningThetag);
260 
261  fHistManager.CreateTH2(Form("hZgPartLevelClosureNoRespFine_%d", cent), Form("Zg response at particle level (closure test, jets not used for the response matrix), %d centrality bin", cent), binEdgesZg.GetSize() - 1, binEdgesZg.GetArray(), binEdgesPtFine.GetSize() - 1, binEdgesPtFine.GetArray());
262  fHistManager.CreateTH2(Form("hZgDetLevelClosureNoRespFine_%d", cent), Form("Zg response at detector level (closure test, jets not used for the response matrix), %d centrality bin", cent), binEdgesZg.GetSize() - 1, binEdgesZg.GetArray(), binEdgesPtFine.GetSize() - 1, binEdgesPtFine.GetArray());
263  fHistManager.CreateTH2(Form("hRgPartLevelClosureNoRespFine_%d", cent), Form("Rg response at particle level (closure test, jets not used for the response matrix), %d centrality bin", cent), binEdgesRg.GetSize() - 1, binEdgesRg.GetArray(), binEdgesPtFine.GetSize() - 1, binEdgesPtFine.GetArray());
264  fHistManager.CreateTH2(Form("hRgDetLevelClosureNoRespFine_%d", cent), Form("Rg response at detector level (closure test, jets not used for the response matrix), %d centrality bin", cent), binEdgesRg.GetSize() - 1, binEdgesRg.GetArray(), binEdgesPtFine.GetSize() - 1, binEdgesPtFine.GetArray());
265  fHistManager.CreateTH2(Form("hNsdPartLevelClosureNoRespFine_%d", cent), Form("Nsd response at particle level (closure test, jets not used for the response matrix), %d centrality bin", cent), binEdgesNsd.GetSize() - 1, binEdgesNsd.GetArray(), binEdgesPtFine.GetSize() - 1, binEdgesPtFine.GetArray());
266  fHistManager.CreateTH2(Form("hNsdDetLevelClosureNoRespFine_%d", cent), Form("Nsd response at detector level (closure test, jets not used for the response matrix), %d centrality bin", cent), binEdgesNsd.GetSize() - 1, binEdgesNsd.GetArray(), binEdgesPtFine.GetSize() - 1, binEdgesPtFine.GetArray());
267  fHistManager.CreateTH2(Form("hThetagPartLevelClosureNoRespFine_%d", cent), Form("Thetag response at particle level (closure test, jets not used for the response matrix), %d centrality bin", cent), binEdgesThetag.GetSize() - 1, binEdgesThetag.GetArray(), binEdgesPtFine.GetSize() - 1, binEdgesPtFine.GetArray());
268  fHistManager.CreateTH2(Form("hThetagDetLevelClosureNoRespFine_%d", cent), Form("Thetag response at detector level (closure test, jets not used for the response matrix), %d centrality bin", cent), binEdgesThetag.GetSize() - 1, binEdgesThetag.GetArray(), binEdgesPtFine.GetSize() - 1, binEdgesPtFine.GetArray());
269  }
270  }
271  }
272  else
273  {
275  fHistManager.CreateTH2("hZgDetLevel", "Zg response at detector level", binEdgesZg.GetSize() - 1, binEdgesZg.GetArray(), binEdgesPtDet.GetSize() - 1, binEdgesPtDet.GetArray());
276  fHistManager.CreateTH2("hZgPartLevel", "Zg response at particle level", binEdgesZg.GetSize() - 1, binEdgesZg.GetArray(), binEdgesPtPart.GetSize() - 1, binEdgesPtPart.GetArray());
277  fHistManager.CreateTH2("hZgPartLevelTruncated", "Zg response at particle level (truncated)", binEdgesZg.GetSize() - 1, binEdgesZg.GetArray(), binEdgesPtPart.GetSize() - 1, binEdgesPtPart.GetArray());
278  fHistManager.CreateTH2("hRgDetLevel", "Rg response at detector level", binEdgesRg.GetSize() - 1, binEdgesRg.GetArray(), binEdgesPtDet.GetSize() - 1, binEdgesPtDet.GetArray());
279  fHistManager.CreateTH2("hRgPartLevel", "Rg response at particle level", binEdgesRg.GetSize() - 1, binEdgesRg.GetArray(), binEdgesPtPart.GetSize() - 1, binEdgesPtPart.GetArray());
280  fHistManager.CreateTH2("hRgPartLevelTruncated", "Rg response at particle level (truncated)", binEdgesRg.GetSize() - 1, binEdgesRg.GetArray(), binEdgesPtPart.GetSize() - 1, binEdgesPtPart.GetArray());
281  fHistManager.CreateTH2("hNsdDetLevel", "Nsd response at detector level", binEdgesNsd.GetSize() - 1, binEdgesNsd.GetArray(), binEdgesPtDet.GetSize() - 1, binEdgesPtDet.GetArray());
282  fHistManager.CreateTH2("hNsdPartLevel", "Nsd response at particle level", binEdgesNsd.GetSize() - 1, binEdgesNsd.GetArray(), binEdgesPtPart.GetSize() - 1, binEdgesPtPart.GetArray());
283  fHistManager.CreateTH2("hNsdPartLevelTruncated", "Nsd response at particle level (truncated)", binEdgesNsd.GetSize() - 1, binEdgesNsd.GetArray(), binEdgesPtPart.GetSize() - 1, binEdgesPtPart.GetArray());
284  fHistManager.CreateTH2("hThetagDetLevel", "Thetag response at detector level", binEdgesThetag.GetSize() - 1, binEdgesThetag.GetArray(), binEdgesPtDet.GetSize() - 1, binEdgesPtDet.GetArray());
285  fHistManager.CreateTH2("hThetagPartLevel", "Thetag response at particle level", binEdgesThetag.GetSize() - 1, binEdgesThetag.GetArray(), binEdgesPtPart.GetSize() - 1, binEdgesPtPart.GetArray());
286  fHistManager.CreateTH2("hThetagPartLevelTruncated", "Thetag response at particle level (truncated)", binEdgesThetag.GetSize() - 1, binEdgesThetag.GetArray(), binEdgesPtPart.GetSize() - 1, binEdgesPtPart.GetArray());
287 
288  // For closure test
289  fHistManager.CreateTH2("hZgPartLevelClosureNoResp", "Zg response at particle level (closure test, jets not used for the response matrix)", binEdgesZg.GetSize() - 1, binEdgesZg.GetArray(), binEdgesPtPart.GetSize() - 1, binEdgesPtPart.GetArray());
290  fHistManager.CreateTH2("hZgDetLevelClosureNoResp", "Zg response at detector level (closure test, jets not used for the response matrix)", binEdgesZg.GetSize() - 1, binEdgesZg.GetArray(), binEdgesPtDet.GetSize() - 1, binEdgesPtDet.GetArray());
291  fHistManager.CreateTH2("hZgPartLevelClosureResp", "Zg response at particle level (closure test, jets used for the response matrix)", binEdgesZg.GetSize() - 1, binEdgesZg.GetArray(), binEdgesPtPart.GetSize() - 1, binEdgesPtPart.GetArray());
292  fHistManager.CreateTH2("hZgDetLevelClosureResp", "Zg response at detector level (closure test, jets used for the response matrix)", binEdgesZg.GetSize() - 1, binEdgesZg.GetArray(), binEdgesPtDet.GetSize() - 1, binEdgesPtDet.GetArray());
293  fHistManager.CreateTH2("hRgPartLevelClosureNoResp", "Rg response at particle level (closure test, jets not used for the response matrix)", binEdgesRg.GetSize() - 1, binEdgesRg.GetArray(), binEdgesPtPart.GetSize() - 1, binEdgesPtPart.GetArray());
294  fHistManager.CreateTH2("hRgDetLevelClosureNoResp", "Rg response at detector level (closure test, jets not used for the response matrix)", binEdgesRg.GetSize() - 1, binEdgesRg.GetArray(), binEdgesPtDet.GetSize() - 1, binEdgesPtDet.GetArray());
295  fHistManager.CreateTH2("hRgPartLevelClosureResp", "Rg response at particle level (closure test, jets used for the response matrix)", binEdgesRg.GetSize() - 1, binEdgesRg.GetArray(), binEdgesPtPart.GetSize() - 1, binEdgesPtPart.GetArray());
296  fHistManager.CreateTH2("hRgDetLevelClosureResp", "Rg response at detector level (closure test, jets used for the response matrix)", binEdgesRg.GetSize() - 1, binEdgesRg.GetArray(), binEdgesPtDet.GetSize() - 1, binEdgesPtDet.GetArray());
297  fHistManager.CreateTH2("hNsdPartLevelClosureNoResp", "Nsd response at particle level (closure test, jets not used for the response matrix)", binEdgesNsd.GetSize() - 1, binEdgesNsd.GetArray(), binEdgesPtPart.GetSize() - 1, binEdgesPtPart.GetArray());
298  fHistManager.CreateTH2("hNsdDetLevelClosureNoResp", "Nsd response at detector level (closure test, jets not used for the response matrix)", binEdgesNsd.GetSize() - 1, binEdgesNsd.GetArray(), binEdgesPtDet.GetSize() - 1, binEdgesPtDet.GetArray());
299  fHistManager.CreateTH2("hNsdPartLevelClosureResp", "Nsd response at particle level (closure test, jets used for the response matrix)", binEdgesNsd.GetSize() - 1, binEdgesNsd.GetArray(), binEdgesPtPart.GetSize() - 1, binEdgesPtPart.GetArray());
300  fHistManager.CreateTH2("hNsdDetLevelClosureResp", "Nsd response at detector level (closure test, jets used for the response matrix)", binEdgesNsd.GetSize() - 1, binEdgesNsd.GetArray(), binEdgesPtDet.GetSize() - 1, binEdgesPtDet.GetArray());
301  fHistManager.CreateTH2("hThetagPartLevelClosureNoResp", "Thetag response at particle level (closure test, jets not used for the response matrix)", binEdgesThetag.GetSize() - 1, binEdgesThetag.GetArray(), binEdgesPtPart.GetSize() - 1, binEdgesPtPart.GetArray());
302  fHistManager.CreateTH2("hThetagDetLevelClosureNoResp", "Thetag response at detector level (closure test, jets not used for the response matrix)", binEdgesThetag.GetSize() - 1, binEdgesThetag.GetArray(), binEdgesPtDet.GetSize() - 1, binEdgesPtDet.GetArray());
303  fHistManager.CreateTH2("hThetagPartLevelClosureResp", "Thetag response at particle level (closure test, jets used for the response matrix)", binEdgesThetag.GetSize() - 1, binEdgesThetag.GetArray(), binEdgesPtPart.GetSize() - 1, binEdgesPtPart.GetArray());
304  fHistManager.CreateTH2("hThetagDetLevelClosureResp", "Thetag response at detector level (closure test, jets used for the response matrix)", binEdgesThetag.GetSize() - 1, binEdgesThetag.GetArray(), binEdgesPtDet.GetSize() - 1, binEdgesPtDet.GetArray());
305 
306  RooUnfoldResponse *r_zg = new RooUnfoldResponse("hZgResponse", "z_{g} response matrix"),
307  *r_rg = new RooUnfoldResponse("hRgResponse", "r_{g} response matrix"),
308  *r_nsd = new RooUnfoldResponse("hNsdResponse", "n_{SD} response matrix"),
309  *r_thetag = new RooUnfoldResponse("hThetagResponse", "#Theta_{g} response matrix");
310  r_zg->Setup((TH1 *)fHistManager.FindObject("hZgDetLevel"), (TH1 *)fHistManager.FindObject("hZgPartLevel"));
311  r_rg->Setup((TH1 *)fHistManager.FindObject("hRgDetLevel"), (TH1 *)fHistManager.FindObject("hRgPartLevel"));
312  r_nsd->Setup((TH1 *)fHistManager.FindObject("hNsdDetLevel"), (TH1 *)fHistManager.FindObject("hNsdPartLevel"));
313  r_thetag->Setup((TH1 *)fHistManager.FindObject("hThetagDetLevel"), (TH1 *)fHistManager.FindObject("hThetagPartLevel"));
314  fZgResponse.push_back(r_zg);
315  fRgResponse.push_back(r_rg);
316  fNsdResponse.push_back(r_nsd);
317  fThetagResponse.push_back(r_thetag);
318  // Create RooUnfold response from THnSparse
319  RooUnfoldResponse *r_zg_closure = new RooUnfoldResponse("hZgResponseClosure", "z_{g} response matrix for the closure test"),
320  *r_rg_closure = new RooUnfoldResponse("hRgResponseClosure", "z_{g} response matrix for the closure test"),
321  *r_nsd_closure = new RooUnfoldResponse("hNsdResponseClosure", "z_{g} response matrix for the closure test"),
322  *r_thetag_closure = new RooUnfoldResponse("hThetagResponseClosure", "#Theta_{g} response matrix for the closure test");
323  r_zg_closure->Setup((TH1 *)fHistManager.FindObject("hZgDetLevel"), (TH1 *)fHistManager.FindObject("hZgPartLevel"));
324  r_rg_closure->Setup((TH1 *)fHistManager.FindObject("hRgDetLevel"), (TH1 *)fHistManager.FindObject("hRgPartLevel"));
325  r_nsd_closure->Setup((TH1 *)fHistManager.FindObject("hNsdDetLevel"), (TH1 *)fHistManager.FindObject("hNsdPartLevel"));
326  r_thetag_closure->Setup((TH1 *)fHistManager.FindObject("hThetagDetLevel"), (TH1 *)fHistManager.FindObject("hThetagPartLevel"));
327  fZgResponseClosure.push_back(r_zg_closure);
328  fRgResponseClosure.push_back(r_rg_closure);
329  fNsdResponseClosure.push_back(r_nsd_closure);
330  fThetagResponseClosure.push_back(r_thetag_closure);
331  }
332 
334  fHistManager.CreateTHnSparse("hZgResponseSparse", "z_{g} response matrix", 4, sparsebinningZg);
335  fHistManager.CreateTHnSparse("hRgResponseSparse", "z_{g} response matrix", 4, sparsebinningRg);
336  fHistManager.CreateTHnSparse("hNsdResponseSparse", "z_{g} response matrix", 4, sparsebinningNsd);
337  fHistManager.CreateTHnSparse("hThetagResponseSparse", "z_{g} response matrix", 4, sparsebinningThetag);
338  fHistManager.CreateTHnSparse("hZgResponseClosureSparse", "z_{g} response matrix for closure test", 4, sparsebinningZg);
339  fHistManager.CreateTHnSparse("hRgResponseClosureSparse", "z_{g} response matrix for closure test", 4, sparsebinningRg);
340  fHistManager.CreateTHnSparse("hNsdResponseClosureSparse", "z_{g} response matrix for closure test", 4, sparsebinningNsd);
341  fHistManager.CreateTHnSparse("hThetagResponseClosureSparse", "z_{g} response matrix for closure test", 4, sparsebinningThetag);
342 
343  fHistManager.CreateTH2("hZgPartLevelClosureNoRespFine", "Zg response at particle level (closure test, jets not used for the response matrix)", binEdgesZg.GetSize() - 1, binEdgesZg.GetArray(), binEdgesPtFine.GetSize() - 1, binEdgesPtFine.GetArray());
344  fHistManager.CreateTH2("hZgDetLevelClosureNoRespFine", "Zg response at detector level (closure test, jets not used for the response matrix)", binEdgesZg.GetSize() - 1, binEdgesZg.GetArray(), binEdgesPtFine.GetSize() - 1, binEdgesPtFine.GetArray());
345  fHistManager.CreateTH2("hRgPartLevelClosureNoRespFine", "Rg response at particle level (closure test, jets not used for the response matrix)", binEdgesRg.GetSize() - 1, binEdgesRg.GetArray(), binEdgesPtFine.GetSize() - 1, binEdgesPtFine.GetArray());
346  fHistManager.CreateTH2("hRgDetLevelClosureNoRespFine", "Rg response at detector level (closure test, jets not used for the response matrix)", binEdgesRg.GetSize() - 1, binEdgesRg.GetArray(), binEdgesPtFine.GetSize() - 1, binEdgesPtFine.GetArray());
347  fHistManager.CreateTH2("hNsdPartLevelClosureNoRespFine", "Nsd response at particle level (closure test, jets not used for the response matrix)", binEdgesNsd.GetSize() - 1, binEdgesNsd.GetArray(), binEdgesPtFine.GetSize() - 1, binEdgesPtFine.GetArray());
348  fHistManager.CreateTH2("hNsdDetLevelClosureNoRespFine", "Nsd response at detector level (closure test, jets not used for the response matrix)", binEdgesNsd.GetSize() - 1, binEdgesNsd.GetArray(), binEdgesPtFine.GetSize() - 1, binEdgesPtFine.GetArray());
349  fHistManager.CreateTH2("hThetagPartLevelClosureNoRespFine", "Thetag response at particle level (closure test, jets not used for the response matrix)", binEdgesThetag.GetSize() - 1, binEdgesThetag.GetArray(), binEdgesPtFine.GetSize() - 1, binEdgesPtFine.GetArray());
350  fHistManager.CreateTH2("hThetagDetLevelClosureNoRespFine", "Thetag response at detector level (closure test, jets not used for the response matrix)", binEdgesThetag.GetSize() - 1, binEdgesThetag.GetArray(), binEdgesPtFine.GetSize() - 1, binEdgesPtFine.GetArray());
351  }
352  }
353 
354  for (auto h : *fHistManager.GetListOfHistograms())
355  fOutput->Add(h);
356  for (auto r : fZgResponse)
357  fOutput->Add(r);
358  for (auto r : fZgResponseClosure)
359  fOutput->Add(r);
360  for (auto r : fRgResponse)
361  fOutput->Add(r);
362  for (auto r : fRgResponseClosure)
363  fOutput->Add(r);
364  for (auto r : fNsdResponse)
365  fOutput->Add(r);
366  for (auto r : fNsdResponseClosure)
367  fOutput->Add(r);
368  for (auto r : fThetagResponse)
369  fOutput->Add(r);
370  for (auto r : fThetagResponseClosure)
371  fOutput->Add(r);
372 
373  PostData(1, fOutput);
374 }
375 
377 {
378  if (!fMCRejectFilter)
379  return true;
380  if (!(fIsPythia || fIsHerwig))
381  return true; // Only relevant for pt-hard production
382  AliDebugStream(1) << "Using custom MC outlier rejection" << std::endl;
384  if (!partjets)
385  return true;
386 
387  // Check whether there is at least one particle level jet with pt above n * event pt-hard
388  auto jetiter = partjets->accepted();
389  auto max = std::max_element(jetiter.begin(), jetiter.end(), [](const AliEmcalJet *lhs, const AliEmcalJet *rhs) { return lhs->Pt() < rhs->Pt(); });
390  if (max != jetiter.end())
391  {
392  // At least one jet found with pt > n * pt-hard
393  AliDebugStream(1) << "Found max jet with pt " << (*max)->Pt() << " GeV/c" << std::endl;
394  if ((*max)->Pt() > fPtHardAndJetPtFactor * fPtHard)
395  return false;
396  }
397  return true;
398 }
399 
401 {
402  enum EPointSD_t {
403  kIndSDDet = 0,
404  kIndPtDet = 1,
405  kIndSDPart = 2,
406  kIndPtPart = 3
407  };
413  double Rjet = detLevelJets->GetJetRadius();
414  if (!(partLevelJets || detLevelJets))
415  {
416  AliErrorStream() << "Either of the jet containers not found" << std::endl;
417  return kFALSE;
418  }
419 
420  if (fSampleFraction < 1.)
421  {
422  if (fSampleTrimmer->Uniform() > fSampleFraction)
423  return false;
424  }
425 
426  // get truncations at detector level
427  TString histname;
428  if (fForceBeamType != kpp)
429  {
430  histname = TString::Format("hZgDetLevel_%d", fCentBin);
431  }
432  else
433  {
434  histname = "hZgDetLevel";
435  }
436  TH2D *fZgDetLevel = (TH2D *)fHistManager.FindObject(histname);
437  auto ptmindet = fHasResponseMatrixRooUnfold ? fZgDetLevel->GetYaxis()->GetBinLowEdge(1) : 0.,
438  ptmaxdet = fHasResponseMatrixRooUnfold ? fZgDetLevel->GetYaxis()->GetBinUpEdge(fZgDetLevel->GetYaxis()->GetNbins()) : 1000.;
439 
440  //when embedding and doing the constituent subtraction there is an additional step to the detector (or hybrid) to particle level matching because the detector jet (or hybrid) is the constituent subtracted jet which is not matched so we need to find the unsubtracted jet that it corresponds to and get the matched jets from there
441  AliJetContainer *jetContUS = nullptr;
442  if (fIsEmbeddedEvent)
444 
445  for (auto detjet : detLevelJets->accepted())
446  {
447  AliEmcalJet *partjet = nullptr;
448  //variables for embedded pbpb data
449  AliEmcalJet *jetUS = nullptr;
450  Int_t ilab = -1;
451  //for embedding, find the unsubtracted jet and get it's matched detector level jet
452  if (fIsEmbeddedEvent)
453  {
454  for (Int_t i = 0; i < jetContUS->GetNJets(); i++)
455  {
456  jetUS = jetContUS->GetJet(i);
457  if (jetUS->GetLabel() == detjet->GetLabel())
458  {
459  ilab = i;
460  break;
461  }
462  }
463  if (ilab == -1)
464  continue;
465  jetUS = jetContUS->GetJet(ilab);
466  partjet = jetUS->ClosestJet();
467  }
468  //if we aren't embedding then just find the matched jet
469  else
470  partjet = detjet->ClosestJet();
471  if (!partjet)
472  continue;
473  //one extra level of matching needed for embedding to go from detector to particle level
474  if (fIsEmbeddedEvent)
475  {
476  partjet = partjet->ClosestJet();
477  if (!partjet)
478  continue;
479  }
480 
481  //cut on the shared pt fraction, when embedding the unsubtracted jet should be used
482  Double_t fraction = 0;
483  if (fIsEmbeddedEvent)
484  fraction = jetContUS->GetFractionSharedPt(jetUS);
485  else
486  fraction = detLevelJets->GetFractionSharedPt(detjet);
487  if (fraction < fMinFractionShared)
488  continue;
489 
490  // sample splitting (for closure test)
491  bool closureUseResponse = (fSampleSplitter->Uniform() < fFractionResponseClosure);
492 
493  // Get the softdrop response
494  std::vector<double> softdropDet, softdropPart;
495 
496  try
497  {
498  softdropDet = MakeSoftdrop(*detjet, detLevelJets->GetJetRadius(), tracks, clusters);
499  softdropPart = MakeSoftdrop(*partjet, partLevelJets->GetJetRadius(), particles, nullptr);
500  bool untaggedDet = softdropDet[0] < fZcut,
501  untaggedPart = softdropPart[0] < fZcut;
502  Double_t pointZg[4] = {softdropDet[0], detjet->Pt(), softdropPart[0], partjet->Pt()},
503  pointRg[4] = {untaggedDet ? -0.01 : softdropDet[2], detjet->Pt(), untaggedPart ? -0.01 : softdropPart[2], partjet->Pt()},
504  pointNsd[4] = {untaggedDet ? -1. : softdropDet[5], detjet->Pt(), untaggedPart ? -1. : softdropPart[5], partjet->Pt()},
505  pointThetag[4] = {untaggedDet ? -0.05 : softdropDet[2]/Rjet, detjet->Pt(), untaggedPart ? -0.05 : softdropPart[2]/Rjet, partjet->Pt()};
506  if (fForceBeamType != kpp)
507  {
509  fHistManager.FillTH1(Form("hZgPartLevel_%d", fCentBin), pointZg[kIndSDPart], pointZg[kIndPtPart]);
510  fHistManager.FillTH1(Form("hRgPartLevel_%d", fCentBin), pointRg[kIndSDPart], pointRg[kIndPtPart]);
511  fHistManager.FillTH1(Form("hNsdPartLevel_%d", fCentBin), pointNsd[kIndSDPart], pointNsd[kIndPtPart]);
512  fHistManager.FillTH1(Form("hThetagPartLevel_%d", fCentBin), pointThetag[kIndSDPart], pointThetag[kIndPtPart]);
513  }
514  }
515  else
516  {
518  fHistManager.FillTH1("hZgPartLevel", pointZg[kIndSDPart], pointZg[kIndPtPart]);
519  fHistManager.FillTH1("hRgPartLevel", pointRg[kIndSDPart], pointRg[kIndPtPart]);
520  fHistManager.FillTH1("hNsdPartLevel", pointNsd[kIndSDPart], pointNsd[kIndPtPart]);
521  fHistManager.FillTH1("hThetagPartLevel", pointThetag[kIndSDPart], pointThetag[kIndPtPart]);
522  }
523  }
524  if (detjet->Pt() >= ptmindet && detjet->Pt() <= ptmaxdet)
525  {
526  if (fForceBeamType != kpp)
527  {
529  fHistManager.FillTH2(Form("hZgPartLevelTruncated_%d", fCentBin), pointZg[kIndSDPart], pointZg[kIndPtPart]);
530  fHistManager.FillTH2(Form("hZgDetLevel_%d", fCentBin), pointZg[kIndSDDet], pointZg[kIndPtDet]);
531  fHistManager.FillTH2(Form("hRgPartLevelTruncated_%d", fCentBin), pointRg[kIndSDPart], pointRg[kIndPtPart]);
532  fHistManager.FillTH2(Form("hRgDetLevel_%d", fCentBin), pointRg[kIndSDDet], pointRg[kIndPtDet]);
533  fHistManager.FillTH2(Form("hNsdPartLevelTruncated_%d", fCentBin), pointNsd[kIndSDPart], pointNsd[kIndPtPart]);
534  fHistManager.FillTH2(Form("hNsdDetLevel_%d", fCentBin), pointNsd[kIndSDDet], pointNsd[kIndPtDet]);
535  fHistManager.FillTH2(Form("hThetagPartLevelTruncated_%d", fCentBin), pointThetag[kIndSDPart], pointThetag[kIndPtPart]);
536  fHistManager.FillTH2(Form("hThetagDetLevel_%d", fCentBin), pointThetag[kIndSDDet], pointThetag[kIndPtDet]);
537  fZgResponse[fCentBin]->Fill(pointZg[kIndSDDet], pointZg[kIndPtDet], pointZg[kIndSDPart], pointZg[kIndPtPart]);
538  fRgResponse[fCentBin]->Fill(pointRg[kIndSDDet], pointRg[kIndPtDet], pointRg[kIndSDPart], pointRg[kIndPtPart]);
539  fNsdResponse[fCentBin]->Fill(pointNsd[kIndSDDet], pointNsd[kIndPtDet], pointNsd[kIndSDPart], pointNsd[kIndPtPart]);
540  fThetagResponse[fCentBin]->Fill(pointThetag[kIndSDDet], pointThetag[kIndPtDet], pointThetag[kIndSDPart], pointThetag[kIndPtPart]);
541  }
543  fHistManager.FillTHnSparse(Form("hZgResponseSparse_%d", fCentBin), pointZg);
544  fHistManager.FillTHnSparse(Form("hRgResponseSparse_%d", fCentBin), pointRg);
545  fHistManager.FillTHnSparse(Form("hNsdResponseSparse_%d", fCentBin), pointNsd);
546  fHistManager.FillTHnSparse(Form("hThetagResponseSparse_%d", fCentBin), pointThetag);
547  }
548  }
549  else
550  {
552  fHistManager.FillTH2("hZgPartLevelTruncated", pointZg[kIndSDPart], pointZg[kIndPtPart]);
553  fHistManager.FillTH2("hZgDetLevel", pointZg[kIndSDDet], pointZg[kIndPtDet]);
554  fHistManager.FillTH2("hRgPartLevelTruncated", pointRg[kIndSDPart], pointRg[kIndPtPart]);
555  fHistManager.FillTH2("hRgDetLevel", pointRg[kIndSDDet], pointRg[kIndPtDet]);
556  fHistManager.FillTH2("hNsdPartLevelTruncated", pointNsd[kIndSDPart], pointNsd[kIndPtPart]);
557  fHistManager.FillTH2("hNsdDetLevel", pointNsd[kIndSDDet], pointNsd[kIndPtDet]);
558  fHistManager.FillTH2("hThetagPartLevelTruncated", pointThetag[kIndSDPart], pointThetag[kIndPtPart]);
559  fHistManager.FillTH2("hThetagDetLevel", pointThetag[kIndSDDet], pointThetag[kIndPtDet]);
560  fZgResponse[0]->Fill(pointZg[kIndSDDet], pointZg[kIndPtDet], pointZg[kIndSDPart], pointZg[kIndPtPart]);
561  fRgResponse[0]->Fill(pointRg[kIndSDDet], pointRg[kIndPtDet], pointRg[kIndSDPart], pointRg[kIndPtPart]);
562  fNsdResponse[0]->Fill(pointNsd[kIndSDDet], pointNsd[kIndPtDet], pointNsd[kIndSDPart], pointNsd[kIndPtPart]);
563  fThetagResponse[0]->Fill(pointThetag[kIndSDDet], pointThetag[kIndPtDet], pointThetag[kIndSDPart], pointThetag[kIndPtPart]);
564  }
566  fHistManager.FillTHnSparse("hZgResponseSparse", pointZg);
567  fHistManager.FillTHnSparse("hRgResponseSparse", pointRg);
568  fHistManager.FillTHnSparse("hNsdResponseSparse", pointNsd);
569  fHistManager.FillTHnSparse("hThetagResponseSparse", pointThetag);
570  }
571  }
572  if (closureUseResponse)
573  {
574  if (fForceBeamType != kpp)
575  {
577  fZgResponseClosure[fCentBin]->Fill(pointZg[kIndSDDet], pointZg[kIndPtDet], pointZg[kIndSDPart], pointZg[kIndPtPart]);
578  fRgResponseClosure[fCentBin]->Fill(pointRg[kIndSDDet], pointRg[kIndPtDet], pointRg[kIndSDPart], pointRg[kIndPtPart]);
579  fNsdResponseClosure[fCentBin]->Fill(pointNsd[kIndSDDet], pointNsd[kIndPtDet], pointNsd[kIndSDPart], pointNsd[kIndPtPart]);
580  fThetagResponseClosure[fCentBin]->Fill(pointThetag[kIndSDDet], pointThetag[kIndPtDet], pointThetag[kIndSDPart], pointThetag[kIndPtPart]);
581  fHistManager.FillTH2(Form("hZgDetLevelClosureResp_%d", fCentBin), pointZg[kIndSDDet], pointZg[kIndPtDet]);
582  fHistManager.FillTH2(Form("hZgPartLevelClosureResp_%d", fCentBin), pointZg[kIndSDPart], pointZg[kIndPtPart]);
583  fHistManager.FillTH2(Form("hRgDetLevelClosureResp_%d", fCentBin), pointRg[kIndSDDet], pointRg[kIndPtDet]);
584  fHistManager.FillTH2(Form("hRgPartLevelClosureResp_%d", fCentBin), pointRg[kIndSDPart], pointRg[kIndPtPart]);
585  fHistManager.FillTH2(Form("hNsdDetLevelClosureResp_%d", fCentBin), pointNsd[kIndSDDet], pointNsd[kIndPtDet]);
586  fHistManager.FillTH2(Form("hNsdPartLevelClosureResp_%d", fCentBin), pointNsd[kIndSDPart], pointNsd[kIndPtPart]);
587  fHistManager.FillTH2(Form("hThetagDetLevelClosureResp_%d", fCentBin), pointThetag[kIndSDDet], pointThetag[kIndPtDet]);
588  fHistManager.FillTH2(Form("hThetagPartLevelClosureResp_%d", fCentBin), pointThetag[kIndSDPart], pointThetag[kIndPtPart]);
589  }
591  fHistManager.FillTHnSparse(Form("hZgResponseClosureSparse_%d", fCentBin), pointZg);
592  fHistManager.FillTHnSparse(Form("hRgResponseClosureSparse_%d", fCentBin), pointRg);
593  fHistManager.FillTHnSparse(Form("hNsdResponseClosureSparse_%d", fCentBin), pointNsd);
594  fHistManager.FillTHnSparse(Form("hThetagResponseClosureSparse_%d", fCentBin), pointThetag);
595  }
596  }
597  else
598  {
600  fZgResponseClosure[0]->Fill(pointZg[kIndSDDet], pointZg[kIndPtDet], pointZg[kIndSDPart], pointZg[kIndPtPart]);
601  fRgResponseClosure[0]->Fill(pointRg[kIndSDDet], pointRg[kIndPtDet], pointRg[kIndSDPart], pointRg[kIndPtPart]);
602  fNsdResponseClosure[0]->Fill(pointNsd[kIndSDDet], pointNsd[kIndPtDet], pointNsd[kIndSDPart], pointNsd[kIndPtPart]);
603  fThetagResponseClosure[0]->Fill(pointThetag[kIndSDDet], pointThetag[kIndPtDet], pointThetag[kIndSDPart], pointThetag[kIndPtPart]);
604  fHistManager.FillTH2("hZgDetLevelClosureResp", pointZg[kIndSDDet], pointZg[kIndPtDet]);
605  fHistManager.FillTH2("hZgPartLevelClosureResp", pointZg[kIndSDPart], pointZg[kIndPtPart]);
606  fHistManager.FillTH2("hRgDetLevelClosureResp", pointRg[kIndSDDet], pointRg[kIndPtDet]);
607  fHistManager.FillTH2("hRgPartLevelClosureResp", pointRg[kIndSDPart], pointRg[kIndPtPart]);
608  fHistManager.FillTH2("hNsdDetLevelClosureResp", pointNsd[kIndSDDet], pointNsd[kIndPtDet]);
609  fHistManager.FillTH2("hNsdPartLevelClosureResp", pointNsd[kIndSDPart], pointNsd[kIndPtPart]);
610  fHistManager.FillTH2("hThetagDetLevelClosureResp", pointThetag[kIndSDDet], pointThetag[kIndPtDet]);
611  fHistManager.FillTH2("hThetagPartLevelClosureResp", pointThetag[kIndSDPart], pointThetag[kIndPtPart]);
612  }
614  fHistManager.FillTHnSparse("hZgResponseClosureSparse", pointZg);
615  fHistManager.FillTHnSparse("hRgResponseClosureSparse", pointRg);
616  fHistManager.FillTHnSparse("hNsdResponseClosureSparse", pointNsd);
617  fHistManager.FillTHnSparse("hThetagResponseClosureSparse", pointThetag);
618  }
619  }
620  }
621  else
622  {
623  if (fForceBeamType != kpp)
624  {
626  fHistManager.FillTH2(Form("hZgPartLevelClosureNoResp_%d", fCentBin), pointZg[kIndSDPart], pointZg[kIndPtPart]);
627  fHistManager.FillTH2(Form("hZgDetLevelClosureNoResp_%d", fCentBin), pointZg[kIndSDDet], pointZg[kIndPtDet]);
628  fHistManager.FillTH2(Form("hRgPartLevelClosureNoResp_%d", fCentBin), pointRg[kIndSDPart], pointRg[kIndPtPart]);
629  fHistManager.FillTH2(Form("hRgDetLevelClosureNoResp_%d", fCentBin), pointRg[kIndSDDet], pointRg[kIndPtDet]);
630  fHistManager.FillTH2(Form("hNsdPartLevelClosureNoResp_%d", fCentBin), pointNsd[kIndSDPart], pointNsd[kIndPtPart]);
631  fHistManager.FillTH2(Form("hNsdDetLevelClosureNoResp_%d", fCentBin), pointNsd[kIndSDDet], pointNsd[kIndPtDet]);
632  fHistManager.FillTH2(Form("hThetagPartLevelClosureNoResp_%d", fCentBin), pointThetag[kIndSDPart], pointThetag[kIndPtPart]);
633  fHistManager.FillTH2(Form("hThetagDetLevelClosureNoResp_%d", fCentBin), pointThetag[kIndSDDet], pointThetag[kIndPtDet]);
634  }
636  fHistManager.FillTH2(Form("hZgPartLevelClosureNoRespFine_%d", fCentBin), pointZg[kIndSDPart], pointZg[kIndPtPart]);
637  fHistManager.FillTH2(Form("hZgDetLevelClosureNoRespFine_%d", fCentBin), pointZg[kIndSDDet], pointZg[kIndPtDet]);
638  fHistManager.FillTH2(Form("hRgPartLevelClosureNoRespFine_%d", fCentBin), pointRg[kIndSDPart], pointRg[kIndPtPart]);
639  fHistManager.FillTH2(Form("hRgDetLevelClosureNoRespFine_%d", fCentBin), pointRg[kIndSDDet], pointRg[kIndPtDet]);
640  fHistManager.FillTH2(Form("hNsdPartLevelClosureNoRespFine_%d", fCentBin), pointNsd[kIndSDPart], pointNsd[kIndPtPart]);
641  fHistManager.FillTH2(Form("hNsdDetLevelClosureNoRespFine_%d", fCentBin), pointNsd[kIndSDDet], pointNsd[kIndPtDet]);
642  fHistManager.FillTH2(Form("hThetagPartLevelClosureNoRespFine_%d", fCentBin), pointThetag[kIndSDPart], pointThetag[kIndPtPart]);
643  fHistManager.FillTH2(Form("hThetagDetLevelClosureNoRespFine_%d", fCentBin), pointThetag[kIndSDDet], pointThetag[kIndPtDet]);
644  }
645  }
646  else
647  {
649  fHistManager.FillTH2("hZgDetLevelClosureNoResp", pointZg[kIndSDDet], pointZg[kIndPtDet]);
650  fHistManager.FillTH2("hZgPartLevelClosureNoResp", pointZg[kIndSDPart], pointZg[kIndPtPart]);
651  fHistManager.FillTH2("hRgDetLevelClosureNoResp", pointRg[kIndSDDet], pointRg[kIndPtDet]);
652  fHistManager.FillTH2("hRgPartLevelClosureNoResp", pointRg[kIndSDPart], pointRg[kIndPtPart]);
653  fHistManager.FillTH2("hNsdDetLevelClosureNoResp", pointNsd[kIndSDDet], pointNsd[kIndPtDet]);
654  fHistManager.FillTH2("hNsdPartLevelClosureNoResp", pointNsd[kIndSDPart], pointNsd[kIndPtPart]);
655  fHistManager.FillTH2("hThetagDetLevelClosureNoResp", pointThetag[kIndSDDet], pointThetag[kIndPtDet]);
656  fHistManager.FillTH2("hThetagPartLevelClosureNoResp", pointThetag[kIndSDPart], pointThetag[kIndPtPart]);
657  }
659  fHistManager.FillTH2("hZgDetLevelClosureNoRespFine", pointZg[kIndSDDet], pointZg[kIndPtDet]);
660  fHistManager.FillTH2("hZgPartLevelClosureNoRespFine", pointZg[kIndSDPart], pointZg[kIndPtPart]);
661  fHistManager.FillTH2("hRgDetLevelClosureNoRespFine", pointRg[kIndSDDet], pointRg[kIndPtDet]);
662  fHistManager.FillTH2("hRgPartLevelClosureNoRespFine", pointRg[kIndSDPart], pointRg[kIndPtPart]);
663  fHistManager.FillTH2("hNsdDetLevelClosureNoRespFine", pointNsd[kIndSDDet], pointNsd[kIndPtDet]);
664  fHistManager.FillTH2("hNsdPartLevelClosureNoRespFine", pointNsd[kIndSDPart], pointNsd[kIndPtPart]);
665  fHistManager.FillTH2("hThetagDetLevelClosureNoRespFine", pointThetag[kIndSDDet], pointThetag[kIndPtDet]);
666  fHistManager.FillTH2("hThetagPartLevelClosureNoRespFine", pointThetag[kIndSDPart], pointThetag[kIndPtPart]);
667  }
668  }
669  }
670  }
671  }
672  catch (...)
673  {
674  AliErrorStream() << "Error in softdrop evaluation - jet will be ignored" << std::endl;
675  continue;
676  }
677  }
678  return kTRUE;
679 }
680 
681 std::vector<double> AliAnalysisTaskEmcalSoftDropResponse::MakeSoftdrop(const AliEmcalJet &jet, double jetradius, const AliParticleContainer *tracks, const AliClusterContainer *clusters) const
682 {
683  const int kClusterOffset = 30000; // In order to handle tracks and clusters in the same index space the cluster index needs and offset, large enough so that there is no overlap with track indices
684  std::vector<fastjet::PseudoJet> constituents;
685  bool isMC = dynamic_cast<const AliMCParticleContainer *>(tracks);
686  AliDebugStream(2) << "Make new jet substrucutre for " << (isMC ? "MC" : "data") << " jet: Number of tracks " << jet.GetNumberOfTracks() << ", clusters " << jet.GetNumberOfClusters() << std::endl;
687  if (tracks && (fUseChargedConstituents || isMC))
688  { // Neutral particles part of particle container in case of MC
689  AliDebugStream(1) << "Jet substructure: Using charged constituents" << std::endl;
690  for (int itrk = 0; itrk < jet.GetNumberOfTracks(); itrk++)
691  {
692  auto track = jet.Track(itrk);
693  if (!track->Charge() && !fUseNeutralConstituents)
694  continue; // Reject neutral constituents in case of using only charged consituents
695  if (track->Charge() && !fUseChargedConstituents)
696  continue; // Reject charged constituents in case of using only neutral consituents
697  fastjet::PseudoJet constituentTrack(track->Px(), track->Py(), track->Pz(), track->E());
698  constituentTrack.set_user_index(jet.TrackAt(itrk));
699  constituents.push_back(constituentTrack);
700  }
701  }
702 
703  if (clusters && fUseNeutralConstituents)
704  {
705  AliDebugStream(1) << "Jet substructure: Using neutral constituents" << std::endl;
706  for (int icl = 0; icl < jet.GetNumberOfClusters(); icl++)
707  {
708  auto cluster = jet.ClusterAt(icl, clusters->GetArray());
709  TLorentzVector clustervec;
710  cluster->GetMomentum(clustervec, fVertex, (AliVCluster::VCluUserDefEnergy_t)clusters->GetDefaultClusterEnergy());
711  fastjet::PseudoJet constituentCluster(clustervec.Px(), clustervec.Py(), clustervec.Pz(), cluster->GetHadCorrEnergy());
712  constituentCluster.set_user_index(jet.ClusterAt(icl) + kClusterOffset);
713  constituents.push_back(constituentCluster);
714  }
715  }
716 
717  AliDebugStream(3) << "Found " << constituents.size() << " constituents for jet with pt=" << jet.Pt() << " GeV/c" << std::endl;
718  if (!constituents.size())
719  {
720  AliErrorStream() << "Jet has 0 constituents." << std::endl;
721  throw 1;
722  }
723  // Redo jet finding on constituents with a
724  fastjet::JetDefinition jetdef(fastjet::antikt_algorithm, jetradius * 2, static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30);
725  fastjet::ClusterSequence jetfinder(constituents, jetdef);
726  std::vector<fastjet::PseudoJet> outputjets = jetfinder.inclusive_jets(0);
727  auto sdjet = outputjets[0];
728  fastjet::contrib::SoftDrop softdropAlgorithm(fBeta, fZcut);
729  softdropAlgorithm.set_verbose_structure(kTRUE);
730  fastjet::JetAlgorithm reclusterizingAlgorithm;
731  switch (fReclusterizer)
732  {
733  case kCAAlgo:
734  reclusterizingAlgorithm = fastjet::cambridge_aachen_algorithm;
735  break;
736  case kKTAlgo:
737  reclusterizingAlgorithm = fastjet::kt_algorithm;
738  break;
739  case kAKTAlgo:
740  reclusterizingAlgorithm = fastjet::antikt_algorithm;
741  break;
742  };
743 #if FASTJET_VERSION_NUMBER >= 30302
744  fastjet::Recluster reclusterizer(reclusterizingAlgorithm, 1, fastjet::Recluster::keep_only_hardest);
745 #else
746  fastjet::contrib::Recluster reclusterizer(reclusterizingAlgorithm, 1, true);
747 #endif
748  softdropAlgorithm.set_reclustering(kTRUE, &reclusterizer);
749  AliDebugStream(4) << "Jet has " << sdjet.constituents().size() << " constituents" << std::endl;
750  auto groomed = softdropAlgorithm(sdjet);
751  auto softdropstruct = groomed.structure_of<fastjet::contrib::SoftDrop>();
752 
753  std::vector<double> result = {softdropstruct.symmetry(),
754  groomed.m(),
755  softdropstruct.delta_R(),
756  groomed.perp(),
757  softdropstruct.mu(),
758  static_cast<double>(softdropstruct.dropped_count())};
759  return result;
760 }
761 
763 {
764  auto binning = new TCustomBinning;
765  binning->SetMinimum(0);
766  switch (fBinningMode)
767  {
768  case kSDModeINT7:
769  {
770  binning->AddStep(20., 20.);
771  binning->AddStep(40., 10.);
772  binning->AddStep(80., 20.);
773  binning->AddStep(120., 40.);
774  binning->AddStep(240., 120.);
775  break;
776  }
777  case kSDModeEJ1:
778  {
779  binning->AddStep(80., 80.);
780  binning->AddStep(140., 10.);
781  binning->AddStep(200., 20.);
782  binning->AddStep(240., 40.);
783  binning->AddStep(400., 160.);
784  break;
785  }
786  case kSDModeEJ2:
787  {
788  binning->AddStep(70., 70.);
789  binning->AddStep(100., 10.);
790  binning->AddStep(140., 20.);
791  binning->AddStep(400., 260.);
792  break;
793  }
794  };
795  return binning;
796 }
797 
799 {
800  auto binning = new TCustomBinning;
801  switch (fBinningMode)
802  {
803  case kSDModeINT7:
804  {
805  binning->SetMinimum(20);
806  binning->AddStep(40., 5.);
807  binning->AddStep(60., 10.);
808  binning->AddStep(80., 20.);
809  binning->AddStep(120., 40.);
810  break;
811  }
812  case kSDModeEJ1:
813  {
814  binning->SetMinimum(80.);
815  binning->AddStep(120., 5.);
816  binning->AddStep(160., 10.);
817  binning->AddStep(200., 20.);
818  break;
819  }
820  case kSDModeEJ2:
821  {
822  binning->SetMinimum(70.);
823  binning->AddStep(100., 5.);
824  binning->AddStep(120., 10.);
825  binning->AddStep(140., 20.);
826  break;
827  }
828  };
829  return binning;
830 }
831 
833 {
834  auto binning = new TCustomBinning;
835  binning->SetMinimum(0.);
836  binning->AddStep(0.1, 0.1);
837  binning->AddStep(0.5, 0.05);
838  return binning;
839 }
840 
842  auto binning = new TCustomBinning;
843  binning->SetMinimum(-0.05); // Negative bins are for untagged jets
844  binning->AddStep(R, 0.05);
845  return binning;
846 }
847 
849 {
850  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
851 
852  Bool_t isAOD(kFALSE);
853  AliInputEventHandler *inputhandler = static_cast<AliInputEventHandler *>(mgr->GetInputEventHandler());
854  if (inputhandler)
855  {
856  if (inputhandler->IsA() == AliAODInputHandler::Class())
857  {
858  std::cout << "Analysing AOD events\n";
859  isAOD = kTRUE;
860  }
861  else
862  {
863  std::cout << "Analysing ESD events\n";
864  }
865  }
866 
867  std::stringstream taskname;
868  taskname << "SoftdropResponsemaker_R" << std::setw(2) << std::setfill('0') << int(jetradius * 10) << trigger;
869  AliAnalysisTaskEmcalSoftDropResponse *responsemaker = new AliAnalysisTaskEmcalSoftDropResponse(taskname.str().data());
870  responsemaker->SelectCollisionCandidates(AliVEvent::kINT7);
871  responsemaker->SetIsEmbeddedEvent(ifembed);
872  mgr->AddTask(responsemaker);
873 
874  TString partcontname(namepartcont);
875  if (partcontname == "usedefault")
876  partcontname = "mcparticles";
877  AliParticleContainer *particles = responsemaker->AddMCParticleContainer(partcontname.Data());
878  //if embedding need to specify that the particles are embedded
879  if (ifembed)
880  particles->SetIsEmbedding(true);
881  particles->SetMinPt(0.);
882  responsemaker->SetNameMCParticleContainer(partcontname.Data());
883 
884  TString partLevelTag("Jet");
885  //if embedding need to specify the tag of the particle level jet from the embedding framework
886  if (ifembed)
887  partLevelTag = "partLevelJets";
888  AliJetContainer *mcjets = responsemaker->AddJetContainer(
889  jettype,
891  recombinationScheme,
892  jetradius,
894  particles, nullptr, partLevelTag);
895  mcjets->SetJetPtCut(0.);
896  mcjets->SetMaxTrackPt(1000.);
897  mcjets->SetName("partjets");
898  responsemaker->SetNamePartLevelJetContainer(mcjets->GetName());
899 
900  AliTrackContainer *tracks(nullptr);
901  if ((jettype == AliJetContainer::kChargedJet) || (jettype == AliJetContainer::kFullJet))
902  {
904  std::cout << "Track container name: " << tracks->GetName() << std::endl;
905  //if embedding need to specify that the tracks are embedded
906  if (ifembed)
907  tracks->SetIsEmbedding(true);
908  tracks->SetMinPt(0.15);
909  }
910  AliClusterContainer *clusters(nullptr);
911  if ((jettype == AliJetContainer::kFullJet) || (jettype == AliJetContainer::kNeutralJet))
912  {
913  std::cout << "Using full or neutral jets ..." << std::endl;
915  std::cout << "Cluster container name: " << clusters->GetName() << std::endl;
916  clusters->SetClusHadCorrEnergyCut(0.3); // 300 MeV E-cut
917  clusters->SetDefaultClusterEnergy(AliVCluster::kHadCorr);
918  }
919  else
920  {
921  std::cout << "Using charged jets ... " << std::endl;
922  }
923 
924  TString detLevelTag("Jet");
925  //if embedding need to specify the tag of the detector (or hybrid) level jet from the embedding framework
926  if (ifembed)
927  detLevelTag = "hybridLevelJets";
928  AliJetContainer *datajets = responsemaker->AddJetContainer(
929  jettype,
931  recombinationScheme,
932  jetradius,
934  tracks, clusters, detLevelTag);
935  datajets->SetJetPtCut(0.);
936  datajets->SetMaxTrackPt(1000.);
937  datajets->SetName("detjets");
938  //if embedding then this jet is the unsubtracted jet and the subtracted jet is added in the run macro, if not then it is the detector level jet
939  if (!ifembed)
940  responsemaker->SetNameDetLevelJetContainer(datajets->GetName());
941  else
942  responsemaker->SetNameUnSubLevelJetContainer(datajets->GetName());
943 
944  std::string jettypestring;
945  switch (jettype)
946  {
948  jettypestring = "FullJets";
949  break;
951  jettypestring = "ChargedJets";
952  break;
954  jettypestring = "NeutralJets";
955  break;
956  default:
957  jettypestring = "Undef";
958  };
959 
960  EBinningMode_t binmode(kSDModeINT7);
961  std::string triggerstring(trigger);
962  if (triggerstring == "EJ1")
963  binmode = kSDModeEJ1;
964  else if (triggerstring == "EJ2")
965  binmode = kSDModeEJ2;
966  responsemaker->SetBinningMode(binmode);
967 
968  // Connecting containers
969  std::stringstream outputfile, histname;
970  outputfile << mgr->GetCommonFileName() << ":SoftDropResponse_" << jettypestring << "_R" << std::setw(2) << std::setfill('0') << int(jetradius * 10.) << "_" << trigger;
971  histname << "SoftDropResponseHistos_" << jettypestring << "_R" << std::setw(2) << std::setfill('0') << int(jetradius * 10.) << "_" << trigger;
972  mgr->ConnectInput(responsemaker, 0, mgr->GetCommonInputContainer());
973  mgr->ConnectOutput(responsemaker, 1, mgr->CreateContainer(histname.str().data(), AliEmcalList::Class(), AliAnalysisManager::kOutputContainer, outputfile.str().data()));
974 
975  return responsemaker;
976 }
Bool_t fHasResponseMatrixSparse
Fill also THnSparse representation of response matrix.
Bool_t fIsPythia
trigger, if it is a PYTHIA production
double Double_t
Definition: External.C:58
std::vector< RooUnfoldResponse * > fNsdResponse
! RooUnfold response for n_sd
Class creating a linear binning, used in the histogram manager.
AliEmcalJet * ClosestJet() const
Definition: AliEmcalJet.h:341
AliJetContainer * GetJetContainer(Int_t i=0) const
void SetName(const char *n)
Set the name of the class of the objets inside the underlying array.
Float_t fMinFractionShared
only fill histos for jets if shared fraction larger than X
std::vector< RooUnfoldResponse * > fThetagResponseClosure
! RooUnfold response for n_sd for the closure test
Container with name, TClonesArray and cuts for particles.
void FillTH2(const char *hname, double x, double y, double weight=1., Option_t *opt="")
Fill a 2D histogram within the container.
Int_t GetLabel() const
Definition: AliEmcalJet.h:124
Int_t ClusterAt(Int_t idx) const
Definition: AliEmcalJet.h:137
void SetMinPt(Double_t min)
AliVParticle * Track(Int_t idx) const
TString fNameDetLevelJetContainer
Name of the detector (or hybrid if embedding) level jet container.
std::vector< RooUnfoldResponse * > fZgResponseClosure
! RooUnfold response for z_g for the closure test
AliJetContainer * AddJetContainer(const char *n, TString defaultCutType, Float_t jetRadius=0.4)
Int_t fCentBin
!event centrality bin
virtual void CreateBinEdges(TArrayD &binedges) const =0
Interface for binnings used by the histogram handler.
Definition: TBinning.h:23
AliClusterContainer * AddClusterContainer(const char *n)
Create new cluster container and attach it to the task.
Container for particles within the EMCAL framework.
Int_t GetDefaultClusterEnergy() const
Double_t GetJetRadius(Int_t i=0) const
Int_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:160
void SetIsEmbedding(Bool_t b)
Set embedding status.
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:139
Double_t fFractionResponseClosure
Fraction of events used for the response matrix in the closure test.
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
TPC fiducial acceptance (each eta edge narrowed by jet R)
Definition: AliEmcalJet.h:68
std::vector< RooUnfoldResponse * > fNsdResponseClosure
! RooUnfold response for n_sd for the closure test
void FillTHnSparse(const char *name, const double *x, double weight=1., Option_t *opt="")
TH2 * CreateTH2(const char *name, const char *title, int nbinsx, double xmin, double xmax, int nbinsy, double ymin, double ymax, Option_t *opt="")
Create a new TH2 within the container.
virtual Bool_t CheckMCOutliers()
Filter the mc tails in pt-hard distributions.
TObject * FindObject(const char *name) const
Find an object inside the container.
int Int_t
Definition: External.C:63
UShort_t GetNumberOfClusters() const
Definition: AliEmcalJet.h:138
void SetJetPtCut(Float_t cut)
TPC acceptance.
Definition: AliEmcalJet.h:67
THashList * GetListOfHistograms() const
Get the list of histograms.
Definition: THistManager.h:671
std::vector< RooUnfoldResponse * > fRgResponse
! RooUnfold response for r_g
Bool_t fIsHerwig
trigger, if it is a HERWIG production
Int_t GetNJets() const
BeamType fForceBeamType
forced beam type
Definition: External.C:228
Int_t fNcentBins
how many centrality bins
AliClusterContainer * GetClusterContainer(Int_t i=0) const
Get cluster container attached to this task.
Helper class creating user defined custom binning.
Float_t fPtHardAndJetPtFactor
Factor between ptHard and jet pT to reject/accept event.
static AliAnalysisTaskEmcalSoftDropResponse * AddTaskEmcalSoftDropResponse(Double_t jetradius, AliJetContainer::EJetType_t jettype, AliJetContainer::ERecoScheme_t recombinationScheme, bool ifembed, const char *namepartcont, const char *trigger)
virtual bool Run()
Run function. This is the core function of the analysis and contains the user code. Therefore users have to implement this function.
TClonesArray * GetArray() const
Bool_t fMCRejectFilter
enable the filtering of events by tail rejection
AliMCParticleContainer * AddMCParticleContainer(const char *n)
Create new container for MC particles and attach it to the task.
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
Fill a 1D histogram within the container.
Double_t Pt() const
Definition: AliEmcalJet.h:109
std::vector< RooUnfoldResponse * > fThetagResponse
! RooUnfold response for theta_g
const char * GetName() const
Bool_t isMC
Float_t fPtHard
!event -hard
Float_t GetJetRadius() const
AliEmcalList * fOutput
!output list
TString fNameUnSubLevelJetContainer
Name of the unsubtracted hybrid level jet container.
Double_t fVertex[3]
!event vertex
AliTrackContainer * AddTrackContainer(const char *n)
Create new track container and attach it to the task.
AliTrackContainer * GetTrackContainer(Int_t i=0) const
void SetMakeGeneralHistograms(Bool_t g)
Enable general histograms.
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:51
std::vector< RooUnfoldResponse * > fRgResponseClosure
! RooUnfold response for r_g for the closure test
Double32_t fSampleFraction
Fraction of statistics used for the analysis.
Double_t GetFractionSharedPt(const AliEmcalJet *jet, AliParticleContainer *cont2=0x0) const
TString fNamePartLevelJetContainer
Name of the particle level jet container.
void UserCreateOutputObjects()
Main initialization function on the worker.
bool Bool_t
Definition: External.C:53
void SetDefaultClusterEnergy(Int_t d)
void SetMaxTrackPt(Float_t b)
THnSparse * CreateTHnSparse(const char *name, const char *title, int ndim, const int *nbins, const double *min, const double *max, Option_t *opt="")
Create a new THnSparse within the container.
std::vector< RooUnfoldResponse * > fZgResponse
! RooUnfold response for z_g
Container structure for EMCAL clusters.
Container for MC-true particles within the EMCAL framework.
EMCal fiducial acceptance (each eta, phi edge narrowed by jet R)
Definition: AliEmcalJet.h:70
Container for jet within the EMCAL jet framework.
Definition: External.C:196
std::vector< double > MakeSoftdrop(const AliEmcalJet &jet, double jetradius, const AliParticleContainer *tracks, const AliClusterContainer *clusters) const
void SetClusHadCorrEnergyCut(Double_t cut)
void SetMinimum(Double_t min)
AliEmcalJet * GetJet(Int_t i) const
static TString TrackContainerNameFactory(Bool_t isAOD)
Get name of the default track container.
static TString ClusterContainerNameFactory(Bool_t isAOD)
Get name of the default cluster container.