AliPhysics  608b256 (608b256)
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 EmcalTriggerJets;
49 
52  fHistos(nullptr),
53  fNameDetectorJets(),
54  fNameParticleJets(),
55  fTriggerSelectionString(),
56  fNameTriggerDecisionContainer("EmcalTriggerDecision"),
57  fFractionResponseClosure(0.8),
58  fSampleSplitter(nullptr)
59 {
60 }
61 
62 AliAnalysisTaskEmcalJetEnergyScale::AliAnalysisTaskEmcalJetEnergyScale(const char *name):
63  AliAnalysisTaskEmcalJet(name, true),
64  fHistos(nullptr),
65  fNameDetectorJets(),
66  fNameParticleJets(),
67  fTriggerSelectionString(),
68  fNameTriggerDecisionContainer("EmcalTriggerDecision"),
69  fFractionResponseClosure(0.8),
70  fSampleSplitter(nullptr)
71 {
72  SetUseAliAnaUtils(true);
73  DefineOutput(1, TList::Class());
74 }
75 
77  if(fHistos) delete fHistos;
79 }
80 
83 
84  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()),
85  subsampleBinning(2, -0.5, 1.5), deltaRbinning(20, 0., 1.), statusbinningEff(3, -0.5, 2.5);
86 
87  const TBinning *diffbinning[3] = {&jetPtBinningPart, &nefbinning, &ptdiffbinning},
88  *corrbinning[6] = {&jetPtBinningPart, &jetPtBinningDet, &nefbinning, &deltaRbinning,&subsampleBinning,&subsampleBinning},
89  *effbinning[3] = {&jetPtBinningPart, &jetPtBinningDet, &statusbinningEff};
90 
91  TCustomBinning jetPtBinningCoarseDet, jetPtBinningCoarsePart;
92  jetPtBinningCoarseDet.SetMinimum(20.);
93  jetPtBinningCoarseDet.AddStep(40., 2.);
94  jetPtBinningCoarseDet.AddStep(60., 5.);
95  jetPtBinningCoarseDet.AddStep(120., 10.);
96  jetPtBinningCoarseDet.AddStep(200., 20.);
97  jetPtBinningCoarsePart.SetMinimum(0);
98  jetPtBinningCoarsePart.AddStep(20., 20.);
99  jetPtBinningCoarsePart.AddStep(80., 10.);
100  jetPtBinningCoarsePart.AddStep(200., 20.);
101  jetPtBinningCoarsePart.AddStep(280., 40.);
102  jetPtBinningCoarsePart.AddStep(500., 220.);
103 
104  fHistos = new THistManager("energyScaleHistos");
105  fHistos->CreateTH1("hEventCounter", "Event counter", 1, 0.5, 1.5);
106  fHistos->CreateTH1("hSpectrumTrueFull", "jet pt spectrum part level, not truncated", jetPtBinningCoarsePart);
107  fHistos->CreateTH1("hSpectrumTrueTruncated", "jet pt spectrum particle level, truncated", jetPtBinningCoarsePart);
108  fHistos->CreateTH2("hJetResponseCoarse", "Response matrix, coarse binning", jetPtBinningCoarseDet, jetPtBinningCoarsePart);
109  fHistos->CreateTHnSparse("hPtDiff", "pt diff det/part", 3, diffbinning, "s");
110  fHistos->CreateTHnSparse("hPtCorr", "Correlation det pt / part pt", 5, corrbinning, "s");
111  fHistos->CreateTHnSparse("hJetfindingEfficiency", "Jet finding efficiency", 3, effbinning, "s");
112  for(auto h : *(fHistos->GetListOfHistograms())) fOutput->Add(h);
113 
114  fSampleSplitter = new TRandom;
115 
116  PostData(1, fOutput);
117 }
118 
120  if(!fMCRejectFilter) return true;
121  if(!(fIsPythia || fIsHerwig)) return true; // Only relevant for pt-hard production
122  AliDebugStream(1) << "Using custom MC outlier rejection" << std::endl;
123  auto partjets = GetJetContainer(fNameParticleJets);
124  if(!partjets) return true;
125 
126  // Check whether there is at least one particle level jet with pt above n * event pt-hard
127  auto jetiter = partjets->accepted();
128  auto max = std::max_element(jetiter.begin(), jetiter.end(), [](const AliEmcalJet *lhs, const AliEmcalJet *rhs ) { return lhs->Pt() < rhs->Pt(); });
129  if(max != jetiter.end()) {
130  // At least one jet found with pt > n * pt-hard
131  AliDebugStream(1) << "Found max jet with pt " << (*max)->Pt() << " GeV/c" << std::endl;
132  if((*max)->Pt() > fPtHardAndJetPtFactor * fPtHard) return false;
133  }
134  return true;
135 }
136 
138  AliDebugStream(1) << "Next event" << std::endl;
139  if(!(fInputHandler->IsEventSelected() & AliVEvent::kINT7)) return false;
141  auto mctrigger = static_cast<PWG::EMCAL::AliEmcalTriggerDecisionContainer *>(fInputEvent->FindListObject(fNameTriggerDecisionContainer));
142  AliDebugStream(1) << "Found trigger decision object: " << (mctrigger ? "yes" : "no") << std::endl;
143  if(!mctrigger){
144  AliErrorStream() << "Trigger decision container with name " << fNameTriggerDecisionContainer << " not found in event - not possible to select EMCAL triggers" << std::endl;
145  return false;
146  }
147  if(!mctrigger->IsEventSelected(fTriggerSelectionString)) return false;
148  }
149  AliDebugStream(1) << "event selected" << std::endl;
150  fHistos->FillTH1("hEventCounter", 1);
151 
152  auto detjets = GetJetContainer(fNameDetectorJets),
154  if(!detjets || !partjets) {
155  AliErrorStream() << "At least one jet container missing, exiting ..." << std::endl;
156  return false;
157  }
158  AliDebugStream(1) << "Have both jet containers: part(" << partjets->GetNAcceptedJets() << "|" << partjets->GetNJets() << "), det(" << detjets->GetNAcceptedJets() << "|" << detjets->GetNJets() << ")" << std::endl;
159 
160  std::vector<AliEmcalJet *> acceptedjets;
161  for(auto detjet : detjets->accepted()){
162  AliDebugStream(2) << "Next jet" << std::endl;
163  acceptedjets.push_back(detjet);
164  auto partjet = detjet->ClosestJet();
165  if(!partjet) {
166  AliDebugStream(2) << "No tagged jet" << std::endl;
167  continue;
168  }
169  Bool_t acceptancematch = false;
170  if (partjet->GetJetAcceptanceType() & detjets->GetAcceptanceType()) acceptancematch = true;
171  TVector3 basevec, tagvec;
172  basevec.SetPtEtaPhi(detjet->Pt(), detjet->Eta(), detjet->Phi());
173  tagvec.SetPtEtaPhi(partjet->Pt(), partjet->Eta(), partjet->Phi());
174  double pointCorr[6] = {partjet->Pt(), detjet->Pt(), detjet->NEF(), basevec.DeltaR(tagvec), acceptancematch ? 1. : 0., fSampleSplitter->Uniform() < fFractionResponseClosure ? 0. : 1.},
175  pointDiff[3] = {partjet->Pt(), detjet->NEF(), (detjet->Pt()-partjet->Pt())/partjet->Pt()};
176  fHistos->FillTHnSparse("hPtDiff", pointDiff);
177  fHistos->FillTHnSparse("hPtCorr", pointCorr);
178  fHistos->FillTH2("hJetResponseCoarse", detjet->Pt(), partjet->Pt());
179  fHistos->FillTH1("hSpectrumTrueFull", partjet->Pt());
180  if(detjet->Pt() >= 20. && detjet->Pt() < 200.) fHistos->FillTH1("hSpectrumTrueTruncated", partjet->Pt());
181  }
182 
183  // efficiency x acceptance: Add histos for all accepted and reconstucted accepted jets
184  for(auto partjet : partjets->accepted()){
185  auto detjet = partjet->ClosestJet();
186  double effvec[3] = {partjet->Pt(), 0., 0.};
187  if(detjet) {
188  // Found a match
189  effvec[1] = detjet->Pt();
190  effvec[2] = 1; // Tagged
191  if(std::find(acceptedjets.begin(), acceptedjets.end(), detjet) != acceptedjets.end()) effvec[2] = 2;
192  }
193  fHistos->FillTHnSparse("hJetfindingEfficiency", effvec);
194  }
195 
196  return true;
197 }
198 
200  const std::array<TString, 8> kEMCALTriggers = {
201  "EJ1", "EJ2", "DJ1", "DJ2", "EG1", "EG2", "DG1", "DG2"
202  };
203  bool isEMCAL = false;
204  for(auto emcaltrg : kEMCALTriggers) {
205  if(triggerstring.Contains(emcaltrg)) {
206  isEMCAL = true;
207  break;
208  }
209  }
210  return isEMCAL;
211 }
212 
214  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
215  if(!mgr){
216  ::Error("EmcalTriggerJets::AliAnalysisTaskEmcalJetEnergyScale::AddTaskJetEnergyScale", "No analysis manager available");
217  return nullptr;
218  }
219 
220  auto inputhandler = mgr->GetInputEventHandler();
221  auto isAOD = inputhandler->IsA() == AliAODInputHandler::Class();
222 
223  std::string jettypename;
225  bool addClusterContainer(false), addTrackContainer(false);
226  switch(jettype){
228  jettypename = "FullJet";
229  acceptance = useDCAL ? AliJetContainer::kDCALfid : AliJetContainer::kEMCALfid;
230  addClusterContainer = addTrackContainer = true;
231  break;
233  jettypename = "ChargedJet";
234  acceptance = AliJetContainer::kTPCfid;
235  addTrackContainer = true;
236  break;
238  jettypename = "NeutralJet";
239  acceptance = useDCAL ? AliJetContainer::kDCALfid : AliJetContainer::kEMCALfid;
240  addClusterContainer = true;
241  break;
243  break;
244  };
245 
246  std::stringstream taskname, tag;
247  tag << jettypename << "_R" << std::setw(2) << std::setfill('0') << int(jetradius * 10.) << "_" << trigger;
248  taskname << "EnergyScaleTask_" << tag.str();
249  AliAnalysisTaskEmcalJetEnergyScale *energyscaletask = new AliAnalysisTaskEmcalJetEnergyScale(taskname.str().data());
250  mgr->AddTask(energyscaletask);
251  energyscaletask->SetTriggerName(trigger);
252 
253  TString partcontname(namepartcont);
254  if(partcontname == "usedefault") partcontname = "mcparticles";
255  auto partcont = energyscaletask->AddMCParticleContainer(partcontname.Data());
256  partcont->SetMinPt(0.);
257 
258  AliClusterContainer *clusters(nullptr);
259  if(addClusterContainer) {
261  switch(jettype){
262  case AliJetContainer::kChargedJet: break; // Silence compiler
264  clusters->SetDefaultClusterEnergy(AliVCluster::kHadCorr);
265  clusters->SetClusHadCorrEnergyCut(0.3);
266  break;
268  clusters->SetDefaultClusterEnergy(AliVCluster::kNonLinCorr);
269  clusters->SetClusNonLinCorrEnergyCut(0.3);
270  break;
272  break;
273  };
274  }
275  AliTrackContainer *tracks(nullptr);
276  if(addTrackContainer) {
278  }
279 
280  auto contpartjet = energyscaletask->AddJetContainer(jettype, AliJetContainer::antikt_algorithm, recoscheme, jetradius,
281  acceptance, partcont, nullptr);
282  contpartjet->SetName("particleLevelJets");
283  energyscaletask->SetNamePartJetContainer("particleLevelJets");
284  std::cout << "Adding particle-level jet container with underling array: " << contpartjet->GetArrayName() << std::endl;
285 
286  auto contdetjet = energyscaletask->AddJetContainer(jettype, AliJetContainer::antikt_algorithm, recoscheme, jetradius,
287  acceptance, tracks, clusters);
288  contdetjet->SetName("detectorLevelJets");
289  energyscaletask->SetNameDetJetContainer("detectorLevelJets");
290  std::cout << "Adding detector-level jet container with underling array: " << contdetjet->GetArrayName() << std::endl;
291 
292  std::stringstream outnamebuilder, listnamebuilder;
293  listnamebuilder << "EnergyScaleHists_" << tag.str();
294  outnamebuilder << mgr->GetCommonFileName() << ":EnergyScaleResults_" << tag.str();
295 
296  mgr->ConnectInput(energyscaletask, 0, mgr->GetCommonInputContainer());
297  mgr->ConnectOutput(energyscaletask, 1, mgr->CreateContainer(listnamebuilder.str().data(), TList::Class(), AliAnalysisManager::kOutputContainer, outnamebuilder.str().data()));
298  return energyscaletask;
299 }
EMCal fiducial acceptance (each eta, phi edge narrowed by jet R)
static AliAnalysisTaskEmcalJetEnergyScale * AddTaskJetEnergyScale(AliJetContainer::EJetType_t jetType, AliJetContainer::ERecoScheme_t recoscheme, Double_t radius, Bool_t useDCAL, const char *namepartcont, const char *trigger)
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)
void AddStep(Double_t max, Double_t binwidth)
Interface for binnings used by the histogram handler.
Definition: TBinning.h:23
TString fNameTriggerDecisionContainer
Global trigger decision container.
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)
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.
Double_t fFractionResponseClosure
Fraction of jets used for response in closure test.
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
void SetClusNonLinCorrEnergyCut(Double_t cut)
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 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.
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.
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.
virtual Bool_t CheckMCOutliers()
Filter the mc tails in pt-hard distributions.
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.
void SetClusHadCorrEnergyCut(Double_t cut)
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.