AliPhysics  09a22ae (09a22ae)
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 
37 #include "AliAODInputHandler.h"
38 #include "AliAnalysisManager.h"
42 #include "AliLog.h"
43 #include "AliVEventHandler.h"
44 
46 
47 using namespace EmcalTriggerJets;
48 
51  fHistos(nullptr),
52  fNameDetectorJets(),
53  fNameParticleJets(),
54  fTriggerSelectionString(),
55  fNameTriggerDecisionContainer("EmcalTriggerDecision")
56 {
57 }
58 
59 AliAnalysisTaskEmcalJetEnergyScale::AliAnalysisTaskEmcalJetEnergyScale(const char *name):
60  AliAnalysisTaskEmcalJet(name, true),
61  fHistos(nullptr),
62  fNameDetectorJets(),
63  fNameParticleJets(),
64  fTriggerSelectionString(),
65  fNameTriggerDecisionContainer("EmcalTriggerDecision")
66 {
67  SetUseAliAnaUtils(true);
68  DefineOutput(1, TList::Class());
69 }
70 
72  if(fHistos) delete fHistos;
73 }
74 
77 
78  TLinearBinning jetPtBinningDet(300, 0., 300.), jetPtBinningPart(500, 0., 500), nefbinning(100, 0., 1.), ptdiffbinning(200, -1., 1.), jetEtaBinning(100, -0.9, 0.9), jetPhiBinning(100, 0., TMath::TwoPi());
79 
80  const TBinning *diffbinning[3] = {&jetPtBinningPart, &nefbinning, &ptdiffbinning},
81  *corrbinning[3] = {&jetPtBinningPart, &jetPtBinningDet, &nefbinning},
82  *effbinning[3] = {&jetPtBinningPart, &jetEtaBinning, &jetPhiBinning};
83 
84  TCustomBinning jetPtBinningCoarseDet, jetPtBinningCoarsePart;
85  jetPtBinningCoarseDet.SetMinimum(20.);
86  jetPtBinningCoarseDet.AddStep(40., 2.);
87  jetPtBinningCoarseDet.AddStep(60., 5.);
88  jetPtBinningCoarseDet.AddStep(120., 10.);
89  jetPtBinningCoarseDet.AddStep(200., 20.);
90  jetPtBinningCoarsePart.SetMinimum(0);
91  jetPtBinningCoarsePart.AddStep(20., 20.);
92  jetPtBinningCoarsePart.AddStep(80., 10.);
93  jetPtBinningCoarsePart.AddStep(200., 20.);
94  jetPtBinningCoarsePart.AddStep(280., 40.);
95  jetPtBinningCoarsePart.AddStep(220., 220.);
96 
97  fHistos = new THistManager("energyScaleHistos");
98  fHistos->CreateTH1("hEventCounter", "Event counter", 1, 0.5, 1.5);
99  fHistos->CreateTH1("hSpectrumTrueFull", "jet pt spectrum part level, not truncated", jetPtBinningCoarsePart);
100  fHistos->CreateTH1("hSpectrumTrueTruncated", "jet pt spectrum particle level, truncated", jetPtBinningCoarsePart);
101  fHistos->CreateTH2("hJetResponseCoarse", "Response matrix, coarse binning", jetPtBinningCoarseDet, jetPtBinningCoarsePart);
102  fHistos->CreateTHnSparse("hPtDiff", "pt diff det/part", 3., diffbinning, "s");
103  fHistos->CreateTHnSparse("hPtCorr", "Correlation det pt / part pt", 3., corrbinning, "s");
104  fHistos->CreateTHnSparse("hPartJetsAccepted", "Accepted particle level jets", 3, effbinning, "s");
105  fHistos->CreateTHnSparse("hPartJetsReconstructed", "Accepted and reconstructed particle level jets", 3, effbinning, "s");
106  for(auto h : *(fHistos->GetListOfHistograms())) fOutput->Add(h);
107 
108  PostData(1, fOutput);
109 }
110 
112  AliDebugStream(1) << "Next event" << std::endl;
113  if(!(fInputHandler->IsEventSelected() & AliVEvent::kINT7)) return false;
115  auto mctrigger = static_cast<PWG::EMCAL::AliEmcalTriggerDecisionContainer *>(fInputEvent->FindListObject(fNameTriggerDecisionContainer));
116  AliDebugStream(1) << "Found trigger decision object: " << (mctrigger ? "yes" : "no") << std::endl;
117  if(!mctrigger){
118  AliErrorStream() << "Trigger decision container with name " << fNameTriggerDecisionContainer << " not found in event - not possible to select EMCAL triggers" << std::endl;
119  return false;
120  }
121  if(!mctrigger->IsEventSelected(fTriggerSelectionString)) return false;
122  }
123  AliDebugStream(1) << "event selected" << std::endl;
124  fHistos->FillTH1("hEventCounter", 1);
125 
126  auto detjets = GetJetContainer(fNameDetectorJets),
128  if(!detjets || !partjets) {
129  AliErrorStream() << "At least one jet container missing, exiting ..." << std::endl;
130  return false;
131  }
132  AliDebugStream(1) << "Have both jet containers: part(" << partjets->GetNAcceptedJets() << "|" << partjets->GetNJets() << "), det(" << detjets->GetNAcceptedJets() << "|" << detjets->GetNJets() << ")" << std::endl;
133 
134  std::vector<AliEmcalJet *> taggedjets;
135  for(auto detjet : detjets->accepted()){
136  AliDebugStream(2) << "Next jet" << std::endl;
137  auto partjet = detjet->ClosestJet();
138  if(!partjet) {
139  AliDebugStream(2) << "No tagged jet" << std::endl;
140  continue;
141  }
142  taggedjets.emplace_back(partjet);
143  double pointCorr[3] = {partjet->Pt(), detjet->Pt(), detjet->NEF()},
144  pointDiff[3] = {partjet->Pt(), detjet->NEF(), (detjet->Pt()-partjet->Pt())/partjet->Pt()};
145  fHistos->FillTHnSparse("hPtDiff", pointDiff);
146  fHistos->FillTHnSparse("hPtCorr", pointCorr);
147  fHistos->FillTH2("hJetResponseCoarse", detjet->Pt(), partjet->Pt());
148  fHistos->FillTH1("hSpectrumTrueFull", partjet->Pt());
149  if(detjet->Pt() >= 20. && detjet->Pt() < 200.) fHistos->FillTH1("hSpectrumTrueTruncated", partjet->Pt());
150  }
151 
152  // efficiency x acceptance: Add histos for all accepted and reconstucted accepted jets
153  for(auto partjet : partjets->accepted()){
154  double pvect[3] = {partjet->Pt(), partjet->Eta(), partjet->Phi()};
155  if(pvect[2] < 0) pvect[2] += TMath::TwoPi();
156  fHistos->FillTHnSparse("hPartJetsAccepted", pvect);
157  if(std::find(taggedjets.begin(), taggedjets.end(), partjet) != taggedjets.end()) fHistos->FillTHnSparse("hPartJetsReconstructed", pvect);
158  }
159  return true;
160 }
161 
163  const std::array<TString, 8> kEMCALTriggers = {
164  "EJ1", "EJ2", "DJ1", "DJ2", "EG1", "EG2", "DG1", "DG2"
165  };
166  bool isEMCAL = false;
167  for(auto emcaltrg : kEMCALTriggers) {
168  if(triggerstring.Contains(emcaltrg)) {
169  isEMCAL = true;
170  break;
171  }
172  }
173  return isEMCAL;
174 }
175 
177  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
178  if(!mgr){
179  ::Error("EmcalTriggerJets::AliAnalysisTaskEmcalJetEnergyScale::AddTaskJetEnergyScale", "No analysis manager available");
180  return nullptr;
181  }
182 
183  auto inputhandler = mgr->GetInputEventHandler();
184  auto isAOD = inputhandler->IsA() == AliAODInputHandler::Class();
185 
186  std::string jettypename;
188  bool addClusterContainer(false), addTrackContainer(false);
189  switch(jettype){
191  jettypename = "FullJet";
192  acceptance = useDCAL ? AliJetContainer::kDCALfid : AliJetContainer::kEMCALfid;
193  addClusterContainer = addTrackContainer = true;
194  break;
196  jettypename = "ChargedJet";
197  acceptance = AliJetContainer::kTPCfid;
198  addTrackContainer = true;
199  break;
201  jettypename = "NeutralJet";
202  acceptance = useDCAL ? AliJetContainer::kDCALfid : AliJetContainer::kEMCALfid;
203  addClusterContainer = true;
204  break;
205  };
206 
207  std::stringstream taskname, tag;
208  tag << jettypename << "_R" << std::setw(2) << std::setfill('0') << int(jetradius * 10.) << "_" << trigger;
209  taskname << "EnergyScaleTask_" << tag.str();
210  AliAnalysisTaskEmcalJetEnergyScale *energyscaletask = new AliAnalysisTaskEmcalJetEnergyScale(taskname.str().data());
211  mgr->AddTask(energyscaletask);
212  energyscaletask->SetTriggerName(trigger);
213 
214  auto partcont = energyscaletask->AddMCParticleContainer("mcparticles");
215  partcont->SetMinPt(0.);
216 
217  AliClusterContainer *clusters(nullptr);
218  if(addClusterContainer) {
220  switch(jettype){
221  case AliJetContainer::kChargedJet: break; // Silence compiler
223  clusters->SetDefaultClusterEnergy(AliVCluster::kHadCorr);
224  clusters->SetClusHadCorrEnergyCut(0.3);
225  break;
227  clusters->SetDefaultClusterEnergy(AliVCluster::kNonLinCorr);
228  clusters->SetClusNonLinCorrEnergyCut(0.3);
229  break;
230  };
231  }
232  AliTrackContainer *tracks(nullptr);
233  if(addTrackContainer) {
235  }
236 
237  auto contpartjet = energyscaletask->AddJetContainer(jettype, AliJetContainer::antikt_algorithm, AliJetContainer::E_scheme, jetradius,
238  AliJetContainer::kTPCfid, partcont, nullptr);
239  contpartjet->SetName("particleLevelJets");
240  energyscaletask->SetNamePartJetContainer("particleLevelJets");
241  std::cout << "Adding particle-level jet container with underling array: " << contpartjet->GetArrayName() << std::endl;
242 
243  auto contdetjet = energyscaletask->AddJetContainer(jettype, AliJetContainer::antikt_algorithm, AliJetContainer::E_scheme, jetradius,
244  acceptance, tracks, clusters);
245  contdetjet->SetName("detectorLevelJets");
246  energyscaletask->SetNameDetJetContainer("detectorLevelJets");
247  std::cout << "Adding detector-level jet container with underling array: " << contdetjet->GetArrayName() << std::endl;
248 
249  std::stringstream outnamebuilder, listnamebuilder;
250  listnamebuilder << "EnergyScaleHists_" << tag.str();
251  outnamebuilder << mgr->GetCommonFileName() << ":EnergyScaleResults_" << tag.str();
252 
253  mgr->ConnectInput(energyscaletask, 0, mgr->GetCommonInputContainer());
254  mgr->ConnectOutput(energyscaletask, 1, mgr->CreateContainer(listnamebuilder.str().data(), TList::Class(), AliAnalysisManager::kOutputContainer, outnamebuilder.str().data()));
255  return energyscaletask;
256 }
EMCal fiducial acceptance (each eta, phi edge narrowed by jet R)
double Double_t
Definition: External.C:58
Class creating a linear binning, used in the histogram manager.
AliJetContainer * GetJetContainer(Int_t i=0) const
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.
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.
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
static AliAnalysisTaskEmcalJetEnergyScale * AddTaskJetEnergyScale(AliJetContainer::EJetType_t jetType, Double_t radius, Bool_t useDCAL, const char *trigger)
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.
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.
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.
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.
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.