AliPhysics  master (3d17d9d)
AliAnalysisTaskEmcalJetEnergyScale.cxx
Go to the documentation of this file.
1 /************************************************************************************
2  * Copyright (C) 2017, 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 <algorithm>
28 #include <array>
29 #include <string>
30 #include <sstream>
31 #include <vector>
32 #include <THistManager.h>
33 #include <TCustomBinning.h>
34 #include <TLinearBinning.h>
35 #include <TCustomBinning.h>
36 #include <TRandom.h>
37 
38 #include "AliAODInputHandler.h"
39 #include "AliAnalysisManager.h"
43 #include "AliLog.h"
44 #include "AliVEventHandler.h"
45 
47 
48 using namespace PWGJE::EMCALJetTasks;
49 
52  fHistos(nullptr),
53  fNameDetectorJets(),
54  fNameParticleJets(),
55  fTriggerSelectionString(),
56  fNameTriggerDecisionContainer("EmcalTriggerDecision"),
57  fFractionResponseClosure(0.8),
58  fFillHSparse(false),
59  fScaleShift(0.),
60  fSampleSplitter(nullptr)
61 {
62 }
63 
64 AliAnalysisTaskEmcalJetEnergyScale::AliAnalysisTaskEmcalJetEnergyScale(const char *name):
65  AliAnalysisTaskEmcalJet(name, true),
66  fHistos(nullptr),
67  fNameDetectorJets(),
68  fNameParticleJets(),
69  fTriggerSelectionString(),
70  fNameTriggerDecisionContainer("EmcalTriggerDecision"),
71  fFractionResponseClosure(0.8),
72  fFillHSparse(false),
73  fScaleShift(0.),
74  fSampleSplitter(nullptr)
75 {
76  SetUseAliAnaUtils(true);
77  DefineOutput(1, TList::Class());
78 }
79 
81  if(fHistos) delete fHistos;
83 }
84 
87 
88  TCustomBinning jetPtBinningCoarseDet, jetPtBinningCoarsePart;
89  jetPtBinningCoarseDet.SetMinimum(20.);
90  jetPtBinningCoarseDet.AddStep(40., 2.);
91  jetPtBinningCoarseDet.AddStep(60., 5.);
92  jetPtBinningCoarseDet.AddStep(120., 10.);
93  jetPtBinningCoarseDet.AddStep(200., 20.);
94  jetPtBinningCoarsePart.SetMinimum(0);
95  jetPtBinningCoarsePart.AddStep(20., 20.);
96  jetPtBinningCoarsePart.AddStep(80., 10.);
97  jetPtBinningCoarsePart.AddStep(200., 20.);
98  jetPtBinningCoarsePart.AddStep(280., 40.);
99  jetPtBinningCoarsePart.AddStep(500., 220.);
100 
101  fHistos = new THistManager("energyScaleHistos");
102  fHistos->CreateTH1("hEventCounter", "Event counter", 1, 0.5, 1.5);
103  fHistos->CreateTH2("hJetEnergyScale", "Jet Energy scale; p_{t,part} (GeV/c); (p_{t,det} - p_{t,part})/p_{t,part}" , 400, 0., 400., 200, -1., 1.);
104  fHistos->CreateTH2("hJetEnergyScaleDet", "Jet Energy scale (det); p_{t,det} (GeV/c); (p_{t,det} - p_{t,part})/p_{t,part}" , 400, 0., 400., 200, -1., 1.);
105  fHistos->CreateTH2("hJetResponseFine", "Response matrix, fine binning", 350, 0., 350., 800, 0., 800.);
106  fHistos->CreateTH2("hJetResponseFineClosure", "Response matrix, fine binning, for closure test", 350, 0., 350., 800, 0., 800.);
107  fHistos->CreateTH2("hJetResponseFineNoClosure", "Response matrix, fine binning, for closure test", 350, 0., 350., 800, 0., 800.);
108  fHistos->CreateTH1("hJetSpectrumPartAll", "Part level jet pt spectrum ", 800, 0., 800.);
109  if(fFillHSparse){
110  TLinearBinning jetPtBinningDet(350, 0., 350.), jetPtBinningPart(600, 0., 600), nefbinning(100, 0., 1.), ptdiffbinning(200, -1., 1.), jetEtaBinning(100, -0.9, 0.9), jetPhiBinning(100, 0., TMath::TwoPi()),
111  subsampleBinning(2, -0.5, 1.5), deltaRbinning(20, 0., 1.), statusbinningEff(3, -0.5, 2.5);
112 
113  const TBinning *diffbinning[3] = {&jetPtBinningPart, &nefbinning, &ptdiffbinning},
114  *corrbinning[6] = {&jetPtBinningPart, &jetPtBinningDet, &nefbinning, &deltaRbinning,&subsampleBinning,&subsampleBinning},
115  *effbinning[3] = {&jetPtBinningPart, &jetPtBinningDet, &statusbinningEff};
116 
117  fHistos->CreateTHnSparse("hPtDiff", "pt diff det/part", 3, diffbinning, "s");
118  fHistos->CreateTHnSparse("hPtCorr", "Correlation det pt / part pt", 6, corrbinning, "s");
119  fHistos->CreateTHnSparse("hJetfindingEfficiency", "Jet finding efficiency", 3, effbinning, "s");
120  }
121  for(auto h : *(fHistos->GetListOfHistograms())) fOutput->Add(h);
122 
123  fSampleSplitter = new TRandom;
124 
125  PostData(1, fOutput);
126 }
127 
129  if(!fMCRejectFilter) return true;
130  if(!(fIsPythia || fIsHerwig)) return true; // Only relevant for pt-hard production
131  AliDebugStream(1) << "Using custom MC outlier rejection" << std::endl;
132  auto partjets = GetJetContainer(fNameParticleJets);
133  if(!partjets) return true;
134 
135  // Check whether there is at least one particle level jet with pt above n * event pt-hard
136  auto jetiter = partjets->accepted();
137  auto max = std::max_element(jetiter.begin(), jetiter.end(), [](const AliEmcalJet *lhs, const AliEmcalJet *rhs ) { return lhs->Pt() < rhs->Pt(); });
138  if(max != jetiter.end()) {
139  // At least one jet found with pt > n * pt-hard
140  AliDebugStream(1) << "Found max jet with pt " << (*max)->Pt() << " GeV/c" << std::endl;
141  if((*max)->Pt() > fPtHardAndJetPtFactor * fPtHard) return false;
142  }
143  return true;
144 }
145 
147  AliDebugStream(1) << "Next event" << std::endl;
148  if(!(fInputHandler->IsEventSelected() & AliVEvent::kINT7)) return false;
150  auto mctrigger = static_cast<PWG::EMCAL::AliEmcalTriggerDecisionContainer *>(fInputEvent->FindListObject(fNameTriggerDecisionContainer));
151  AliDebugStream(1) << "Found trigger decision object: " << (mctrigger ? "yes" : "no") << std::endl;
152  if(!mctrigger){
153  AliErrorStream() << "Trigger decision container with name " << fNameTriggerDecisionContainer << " not found in event - not possible to select EMCAL triggers" << std::endl;
154  return false;
155  }
156  if(!mctrigger->IsEventSelected(fTriggerSelectionString)) return false;
157  }
158  AliDebugStream(1) << "event selected" << std::endl;
159  fHistos->FillTH1("hEventCounter", 1);
160 
161  auto detjets = GetJetContainer(fNameDetectorJets),
163  if(!detjets || !partjets) {
164  AliErrorStream() << "At least one jet container missing, exiting ..." << std::endl;
165  return false;
166  }
167  AliDebugStream(1) << "Have both jet containers: part(" << partjets->GetNAcceptedJets() << "|" << partjets->GetNJets() << "), det(" << detjets->GetNAcceptedJets() << "|" << detjets->GetNJets() << ")" << std::endl;
168 
169  std::vector<AliEmcalJet *> acceptedjets;
170  for(auto detjet : detjets->accepted()){
171  AliDebugStream(2) << "Next jet" << std::endl;
172  acceptedjets.push_back(detjet);
173  auto partjet = detjet->ClosestJet();
174  if(!partjet) {
175  AliDebugStream(2) << "No tagged jet" << std::endl;
176  continue;
177  }
178  bool isClosure = fSampleSplitter->Uniform() < fFractionResponseClosure;
179  Double_t detpt = detjet->Pt();
180  if(TMath::Abs(fScaleShift) > DBL_EPSILON){
181  detpt += fScaleShift * detpt;
182  }
183  if(fFillHSparse) {
184  Bool_t acceptancematch = false;
185  if (partjet->GetJetAcceptanceType() & detjets->GetAcceptanceType()) acceptancematch = true;
186  TVector3 basevec, tagvec;
187  basevec.SetPtEtaPhi(detjet->Pt(), detjet->Eta(), detjet->Phi());
188  tagvec.SetPtEtaPhi(partjet->Pt(), partjet->Eta(), partjet->Phi());
189  double pointCorr[6] = {partjet->Pt(), detpt, detjet->NEF(), basevec.DeltaR(tagvec), acceptancematch ? 1. : 0., isClosure ? 0. : 1.},
190  pointDiff[3] = {partjet->Pt(), detjet->NEF(), (detpt-partjet->Pt())/partjet->Pt()};
191  fHistos->FillTHnSparse("hPtDiff", pointDiff);
192  fHistos->FillTHnSparse("hPtCorr", pointCorr);
193  }
194  fHistos->FillTH2("hJetResponseFine", detpt, partjet->Pt());
195  fHistos->FillTH1("hJetEnergyScale", partjet->Pt(), (detpt - partjet->Pt())/partjet->Pt());
196  fHistos->FillTH1("hJetEnergyScaleDet", detpt, (detpt - partjet->Pt())/partjet->Pt());
197  // splitting for closure test
198  if(isClosure) {
199  fHistos->FillTH2("hJetResponseFineClosure", detpt, partjet->Pt());
200  } else {
201  fHistos->FillTH2("hJetResponseFineNoClosure", detpt, partjet->Pt());
202  }
203  }
204 
205  // efficiency x acceptance: Add histos for all accepted and reconstucted accepted jets
206  for(auto partjet : partjets->accepted()){
207  if(fFillHSparse){
208  auto detjet = partjet->ClosestJet();
209  double effvec[3] = {partjet->Pt(), 0., 0.};
210  if(detjet) {
211  // Found a match
212  effvec[1] = detjet->Pt();
213  if(TMath::Abs(fScaleShift) > DBL_EPSILON){
214  effvec[1] += fScaleShift * effvec[1];
215  }
216  effvec[2] = 1; // Tagged
217  if(std::find(acceptedjets.begin(), acceptedjets.end(), detjet) != acceptedjets.end()) effvec[2] = 2;
218  }
219  fHistos->FillTHnSparse("hJetfindingEfficiency", effvec);
220  }
221  fHistos->FillTH1("hJetSpectrumPartAll", partjet->Pt());
222  }
223 
224  return true;
225 }
226 
228  const std::array<TString, 8> kEMCALTriggers = {
229  "EJ1", "EJ2", "DJ1", "DJ2", "EG1", "EG2", "DG1", "DG2"
230  };
231  bool isEMCAL = false;
232  for(auto emcaltrg : kEMCALTriggers) {
233  if(triggerstring.Contains(emcaltrg)) {
234  isEMCAL = true;
235  break;
236  }
237  }
238  return isEMCAL;
239 }
240 
241 AliAnalysisTaskEmcalJetEnergyScale *AliAnalysisTaskEmcalJetEnergyScale::AddTaskJetEnergyScale(AliJetContainer::EJetType_t jettype, AliJetContainer::ERecoScheme_t recoscheme, AliVCluster::VCluUserDefEnergy_t energydef, Double_t jetradius, Bool_t useDCAL, const char *namepartcont, const char *trigger, const char *suffix) {
242  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
243  if(!mgr){
244  ::Error("EmcalTriggerJets::AliAnalysisTaskEmcalJetEnergyScale::AddTaskJetEnergyScale", "No analysis manager available");
245  return nullptr;
246  }
247 
248  auto inputhandler = mgr->GetInputEventHandler();
249  auto isAOD = inputhandler->IsA() == AliAODInputHandler::Class();
250 
251  std::string jettypename;
253  AliJetContainer::EJetType_t mcjettype(jettype);
254  bool addClusterContainer(false), addTrackContainer(false);
255  switch(jettype){
257  jettypename = "FullJet";
258  acceptance = useDCAL ? AliJetContainer::kDCALfid : AliJetContainer::kEMCALfid;
259  addClusterContainer = addTrackContainer = true;
260  break;
262  jettypename = "ChargedJet";
263  acceptance = AliJetContainer::kTPCfid;
264  addTrackContainer = true;
265  mcjettype = AliJetContainer::kChargedJet;
266  break;
268  jettypename = "NeutralJet";
269  acceptance = useDCAL ? AliJetContainer::kDCALfid : AliJetContainer::kEMCALfid;
270  addClusterContainer = true;
271  mcjettype = AliJetContainer::kFullJet; // Correct back neutral detector-level jets to full particle level jets
272  break;
274  break;
275  };
276 
277  std::stringstream taskname, tag;
278  tag << jettypename << "_R" << std::setw(2) << std::setfill('0') << int(jetradius * 10.) << "_" << trigger;
279  if(strlen(suffix)) tag << "_" << suffix;
280  taskname << "EnergyScaleTask_" << tag.str();
281  AliAnalysisTaskEmcalJetEnergyScale *energyscaletask = new AliAnalysisTaskEmcalJetEnergyScale(taskname.str().data());
282  mgr->AddTask(energyscaletask);
283  energyscaletask->SetTriggerName(trigger);
284 
285  TString partcontname(namepartcont);
286  if(partcontname == "usedefault") partcontname = "mcparticles";
287  auto partcont = energyscaletask->AddMCParticleContainer(partcontname.Data());
288  partcont->SetMinPt(0.);
289 
290  AliClusterContainer *clusters(nullptr);
291  if(addClusterContainer) {
293  clusters->SetDefaultClusterEnergy(energydef);
294  clusters->SetClusUserDefEnergyCut(energydef, 0.3);
295  }
296  AliTrackContainer *tracks(nullptr);
297  if(addTrackContainer) {
299  }
300 
301  auto contpartjet = energyscaletask->AddJetContainer(mcjettype, AliJetContainer::antikt_algorithm, recoscheme, jetradius,
302  acceptance, partcont, nullptr);
303  contpartjet->SetName("particleLevelJets");
304  energyscaletask->SetNamePartJetContainer("particleLevelJets");
305  std::cout << "Adding particle-level jet container with underling array: " << contpartjet->GetArrayName() << std::endl;
306 
307  auto contdetjet = energyscaletask->AddJetContainer(jettype, AliJetContainer::antikt_algorithm, recoscheme, jetradius,
308  acceptance, tracks, clusters);
309  contdetjet->SetName("detectorLevelJets");
310  energyscaletask->SetNameDetJetContainer("detectorLevelJets");
311  std::cout << "Adding detector-level jet container with underling array: " << contdetjet->GetArrayName() << std::endl;
312 
313  std::stringstream outnamebuilder, listnamebuilder;
314  listnamebuilder << "EnergyScaleHists_" << tag.str();
315  outnamebuilder << mgr->GetCommonFileName() << ":EnergyScaleResults_" << tag.str();
316 
317  mgr->ConnectInput(energyscaletask, 0, mgr->GetCommonInputContainer());
318  mgr->ConnectOutput(energyscaletask, 1, mgr->CreateContainer(listnamebuilder.str().data(), TList::Class(), AliAnalysisManager::kOutputContainer, outnamebuilder.str().data()));
319  return energyscaletask;
320 }
EMCal fiducial acceptance (each eta, phi edge narrowed by jet R)
Bool_t fIsPythia
trigger, if it is a PYTHIA production
double Double_t
Definition: External.C:58
Class creating a linear binning, used in the histogram manager.
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.
Container with name, TClonesArray and cuts for particles.
void SetUseAliAnaUtils(Bool_t b, Bool_t bRejPilup=kTRUE)
void FillTH2(const char *hname, double x, double y, double weight=1., Option_t *opt="")
Fill a 2D histogram within the container.
void SetMinPt(Double_t min)
AliJetContainer * AddJetContainer(const char *n, TString defaultCutType, Float_t jetRadius=0.4)
Double_t fFractionResponseClosure
Fraction of jets used for response in closure test.
void AddStep(Double_t max, Double_t binwidth)
Interface for binnings used by the histogram handler.
Definition: TBinning.h:23
void SetClusUserDefEnergyCut(Int_t t, Double_t cut)
AliClusterContainer * AddClusterContainer(const char *n)
Create new cluster container and attach it to the task.
DCal fiducial acceptance (each eta, phi edge narrowed by jet R)
static AliAnalysisTaskEmcalJetEnergyScale * AddTaskJetEnergyScale(AliJetContainer::EJetType_t jetType, AliJetContainer::ERecoScheme_t recoscheme, AliVCluster::VCluUserDefEnergy_t energydef, Double_t radius, Bool_t useDCAL, const char *namepartcont, const char *trigger, const char *suffix)
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.
JetAcceptanceType
Bit definition for jet geometry acceptance. Defined here for backwards compatibility. This will be removed. Please use AliEmcalJet::JetAcceptanceType in your code.
THashList * GetListOfHistograms() const
Get the list of histograms.
Definition: THistManager.h:671
Bool_t fIsHerwig
trigger, if it is a HERWIG production
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
Create a new TH1 within the container.
Helper class creating user defined custom binning.
Float_t fPtHardAndJetPtFactor
Factor between ptHard and jet pT to reject/accept event.
virtual Bool_t Run()
Run function. This is the core function of the analysis and contains the user code. Therefore users have to implement this function.
Bool_t fMCRejectFilter
enable the filtering of events by tail rejection
virtual Bool_t CheckMCOutliers()
Filter the mc tails in pt-hard distributions.
AliMCParticleContainer * AddMCParticleContainer(const char *n)
Create new container for MC particles and attach it to the task.
TPC fiducial acceptance (each eta edge narrowed by jet R)
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
Float_t fPtHard
!event -hard
AliEmcalList * fOutput
!output list
AliTrackContainer * AddTrackContainer(const char *n)
Create new track container and attach it to the task.
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:51
Container class for histograms.
Definition: THistManager.h:99
void UserCreateOutputObjects()
Main initialization function on the worker.
bool Bool_t
Definition: External.C:53
void SetDefaultClusterEnergy(Int_t d)
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.
Container structure for EMCAL clusters.
TString fNameTriggerDecisionContainer
Global trigger decision container.
void SetMinimum(Double_t min)
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.