AliPhysics  b76e98e (b76e98e)
AliAnalysisTaskEmcalJetEnergySpectrum.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 <sstream>
30 #include <string>
31 #include <vector>
32 
33 #include <THistManager.h>
34 #include <TLinearBinning.h>
35 #include <TVariableBinning.h>
36 
37 #include "AliAODInputHandler.h"
38 #include "AliAnalysisManager.h"
41 #include "AliEmcalJet.h"
44 #include "AliInputEventHandler.h"
45 #include "AliJetContainer.h"
46 #include "AliLog.h"
47 #include "AliVEvent.h"
48 
50 
51 using namespace EmcalTriggerJets;
52 
55  fHistos(nullptr),
56  fIsMC(false),
57  fTriggerSelectionBits(AliVEvent::kAny),
58  fTriggerSelectionString(""),
59  fRequireSubsetMB(false),
60  fMinBiasTrigger(AliVEvent::kAny),
61  fNameTriggerDecisionContainer("EmcalTriggerDecision"),
62  fUseTriggerSelectionForData(false),
63  fUseDownscaleWeight(false),
64  fNameJetContainer("datajets")
65 {
66  SetUseAliAnaUtils(true);
67 }
68 
70  AliAnalysisTaskEmcalJet(name, true),
72  fIsMC(false),
73  fTriggerSelectionBits(AliVEvent::kAny),
75  fRequireSubsetMB(false),
76  fMinBiasTrigger(AliVEvent::kAny),
77  fNameTriggerDecisionContainer("EmcalTriggerDecision"),
79  fUseDownscaleWeight(false),
80  fNameJetContainer("datajets")
81 {
82  SetUseAliAnaUtils(true);
83 }
84 
86  if(fHistos) delete fHistos;
87 }
88 
91 
92  if(!fUserPtBinning.GetSize()) {
93  // binning not set. apply default binning
94  AliInfoStream() << "Using default pt binning";
95  fUserPtBinning.Set(201);
96  double current(0.);
97  for(int istep = 0; istep < 201; istep++) {
98  fUserPtBinning[istep] = current;
99  current += 1;
100  }
101  }
102 
103  TLinearBinning etabinning(100, -1., 1.), phibinning(100., 0., 7.), nefbinning(100, 0., 1.), trgclusterbinning(kTrgClusterN + 1, -0.5, kTrgClusterN -0.5);
104  TVariableBinning jetptbinning(fUserPtBinning);
105  const TBinning *binnings[5] = {&jetptbinning, &etabinning, &phibinning, &nefbinning, &trgclusterbinning};
106  fHistos = new THistManager(Form("Histos_%s", GetName()));
107  fHistos->CreateTH1("hEventCounter", "Event counter histogram", 1, 0.5, 1.5);
108  fHistos->CreateTH1("hClusterCounter", "Event counter histogram", kTrgClusterN, -0.5, kTrgClusterN - 0.5);
109  fHistos->CreateTHnSparse("hJetTHnSparse", "jet thnsparse", 5, binnings, "s");
110  fHistos->CreateTHnSparse("hMaxJetTHnSparse", "jet thnsparse", 5, binnings, "s");
111 
112  for(auto h : *fHistos->GetListOfHistograms()) fOutput->Add(h);
113  PostData(1, fOutput);
114 }
115 
117  auto datajets = this->GetJetContainer(fNameJetContainer);
118  if(!datajets) {
119  AliErrorStream() << "Jet container " << fNameJetContainer << " not found" << std::endl;
120  return false;
121  }
122 
123  auto trgclusters = GetTriggerClusterIndices(fInputEvent->GetFiredTriggerClasses());
124  fHistos->FillTH1("hEventCounter", 1);
125  AliEmcalJet *maxjet(nullptr);
126  for(auto t : trgclusters) fHistos->FillTH1("hClusterCounter", t);
127  for(auto j : datajets->accepted()){
128  if(!maxjet || (j->E() > maxjet->E())) maxjet = j;
129  double datapoint[5] = {j->Pt(), j->Eta(), j->Phi(), j->NEF(), 0.};
130  for(auto t : trgclusters){
131  datapoint[4] = static_cast<double>(t);
132  fHistos->FillTHnSparse("hJetTHnSparse", datapoint);
133  }
134  }
135 
136  double maxdata[5];
137  memset(maxdata, 0., sizeof(double) * 5);
138  if(maxjet){
139  maxdata[0] = maxjet->Pt();
140  maxdata[1] = maxjet->Eta();
141  maxdata[2] = maxjet->Phi();
142  maxdata[3] = maxjet->NEF();
143  for(auto t : trgclusters){
144  maxdata[4] = static_cast<double>(t);
145  fHistos->FillTHnSparse("hMaxJetTHnSparse", maxdata);
146  }
147  }
148  return true;
149 }
150 
151 std::vector<AliAnalysisTaskEmcalJetEnergySpectrum::TriggerCluster_t> AliAnalysisTaskEmcalJetEnergySpectrum::GetTriggerClusterIndices(const TString &triggerstring) const {
152  // decode trigger string in order to determine the trigger clusters
153  std::vector<TriggerCluster_t> result;
154  result.emplace_back(kTrgClusterANY); // cluster ANY always included
155  if(!fIsMC){
156  // Data - separate trigger clusters
157  std::vector<std::string> clusternames;
158  auto triggerinfos = PWG::EMCAL::Triggerinfo::DecodeTriggerString(triggerstring.Data());
159  for(auto t : triggerinfos) {
160  if(std::find(clusternames.begin(), clusternames.end(), t.Triggercluster()) == clusternames.end()) clusternames.emplace_back(t.Triggercluster());
161  }
162  bool isCENT = (std::find(clusternames.begin(), clusternames.end(), "CENT") != clusternames.end()),
163  isCENTNOTRD = (std::find(clusternames.begin(), clusternames.end(), "CENTNOTRD") != clusternames.end()),
164  isCALO = (std::find(clusternames.begin(), clusternames.end(), "CALO") != clusternames.end()),
165  isCALOFAST = (std::find(clusternames.begin(), clusternames.end(), "CALOFAST") != clusternames.end());
166  if(isCENT || isCENTNOTRD) {
167  if(isCENT) {
168  result.emplace_back(kTrgClusterCENT);
169  if(isCENTNOTRD) {
170  result.emplace_back(kTrgClusterCENTNOTRD);
171  result.emplace_back(kTrgClusterCENTBOTH);
172  } else result.emplace_back(kTrgClusterOnlyCENT);
173  } else {
174  result.emplace_back(kTrgClusterCENTNOTRD);
175  result.emplace_back(kTrgClusterOnlyCENTNOTRD);
176  }
177  }
178  if(isCALO || isCALOFAST) {
179  if(isCALO) {
180  result.emplace_back(kTrgClusterCALO);
181  if(isCALOFAST) {
182  result.emplace_back(kTrgClusterCALOFAST);
183  result.emplace_back(kTrgClusterCALOBOTH);
184  } else result.emplace_back(kTrgClusterOnlyCALO);
185  } else {
186  result.emplace_back(kTrgClusterCALOFAST);
187  result.emplace_back(kTrgClusterOnlyCALOFAST);
188  }
189  }
190  }
191  return result;
192 }
193 
195  if(!fIsMC){
196  // Pure data - do EMCAL trigger selection from selection string
197  if(!(fInputHandler->IsEventSelected() & fTriggerSelectionBits)) return false;
198  if(fTriggerSelectionString.Length()) {
199  if(!fInputEvent->GetFiredTriggerClasses().Contains(fTriggerSelectionString)) return false;
200  if(fRequireSubsetMB && !(fInputHandler->IsEventSelected() & fMinBiasTrigger)) return false; // Require EMCAL trigger to be subset of the min. bias trigger (for efficiency studies)
201  if((fTriggerSelectionString.Contains("EJ") || fTriggerSelectionString.Contains("EG") || fTriggerSelectionString.Contains("DJ") || fTriggerSelectionString.Contains("DG")) && fUseTriggerSelectionForData) {
202  auto trgselresult = static_cast<PWG::EMCAL::AliEmcalTriggerDecisionContainer *>(fInputEvent->FindListObject(fNameTriggerDecisionContainer));
203  AliDebugStream(1) << "Found trigger decision object: " << (trgselresult ? "yes" : "no") << std::endl;
204  if(!trgselresult){
205  AliErrorStream() << "Trigger decision container with name " << fNameTriggerDecisionContainer << " not found in event - not possible to select EMCAL triggers" << std::endl;
206  return false;
207  }
208  if(!trgselresult->IsEventSelected(fTriggerSelectionString)) return false;
209  }
210  }
211  } else {
212  if(!(fInputHandler->IsEventSelected() & AliVEvent::kINT7)) return false;
214  // Simulation - do EMCAL trigger selection from trigger selection object
215  auto mctrigger = static_cast<PWG::EMCAL::AliEmcalTriggerDecisionContainer *>(fInputEvent->FindListObject(fNameTriggerDecisionContainer));
216  AliDebugStream(1) << "Found trigger decision object: " << (mctrigger ? "yes" : "no") << std::endl;
217  if(!mctrigger){
218  AliErrorStream() << "Trigger decision container with name " << fNameTriggerDecisionContainer << " not found in event - not possible to select EMCAL triggers" << std::endl;
219  return false;
220  }
221  if(!mctrigger->IsEventSelected(fTriggerSelectionString)) return false;
222  }
223  }
224  return true;
225 }
226 
227 bool AliAnalysisTaskEmcalJetEnergySpectrum::IsSelectEmcalTriggers(const std::string &triggerstring) const {
228  const std::array<std::string, 8> kEMCALTriggers = {
229  "EJ1", "EJ2", "DJ1", "DJ2", "EG1", "EG2", "DG1", "DG2"
230  };
231  bool isEMCAL = false;
232  for(auto emcaltrg : kEMCALTriggers) {
233  if(triggerstring.find(emcaltrg) != std::string::npos) {
234  isEMCAL = true;
235  break;
236  }
237  }
238  return isEMCAL;
239 }
240 
241 
243  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
244  if(!mgr) {
245  std::cerr << "Analysis manager not initialized" << std::endl;
246  return nullptr;
247  }
248 
249  Bool_t isAOD(kFALSE);
250  AliInputEventHandler *inputhandler = static_cast<AliInputEventHandler *>(mgr->GetInputEventHandler());
251  if(inputhandler) {
252  if(inputhandler->IsA() == AliAODInputHandler::Class()){
253  std::cout << "Analysing AOD events\n";
254  isAOD = kTRUE;
255  } else {
256  std::cout << "Analysing ESD events\n";
257  }
258  }
259 
260  std::string jettypestring;
262  switch(jettype){
263  case AliJetContainer::kChargedJet: jettypestring = "ChargedJets"; acctype = AliJetContainer::kTPCfid; break;
264  case AliJetContainer::kFullJet: jettypestring = "FullJets"; acctype = AliJetContainer::kEMCALfid; break;
265  case AliJetContainer::kNeutralJet: jettypestring = "NeutralJets"; acctype = AliJetContainer::kEMCALfid; break;
266  };
267 
268  std::stringstream tag, outfilename;
269  tag << jettypestring << "_R" << std::setw(2) << std::setfill('0') << int(radius * 10.) << "_" << trigger;
270  if(strlen(suffix)) {
271  tag << "_" << suffix;
272  }
273  auto task = new AliAnalysisTaskEmcalJetEnergySpectrum(Form("JetEnergySpectrum_%s", tag.str().data()));
274  task->SetIsMC(isMC);
275  mgr->AddTask(task);
276 
277  auto contains = [](const std::string &str, const std::string &test) {
278  return str.find(test) != std::string::npos;
279  };
280 
281  std::string trgstr(trigger);
282  if(contains(trgstr, "INT7")) task->SetTriggerSelection(AliVEvent::kINT7, "INT7");
283  else if(contains(trgstr, "EJ1")) task->SetTriggerSelection(AliVEvent::kEMCEJE, "EJ1");
284  else if(contains(trgstr, "EJ2")) task->SetTriggerSelection(AliVEvent::kEMCEJE, "EJ2");
285  else if(contains(trgstr, "EG1")) task->SetTriggerSelection(AliVEvent::kEMCEGA, "EG1");
286  else if(contains(trgstr, "EG2")) task->SetTriggerSelection(AliVEvent::kEMCEGA, "EG2");
287 
288  // Connect particle and cluster container
289  AliTrackContainer *tracks(nullptr);
290  AliClusterContainer *clusters(nullptr);
291  if(jettype == AliJetContainer::kChargedJet || jettype == AliJetContainer::kFullJet) {
293  tracks->SetMinPt(0.15);
294  }
295  if(jettype == AliJetContainer::kNeutralJet || jettype == AliJetContainer::kFullJet){
296  clusters = task->AddClusterContainer(EMCalTriggerPtAnalysis::AliEmcalAnalysisFactory::ClusterContainerNameFactory(isAOD));
297  clusters->SetDefaultClusterEnergy(AliVCluster::kHadCorr);
298  clusters->SetClusHadCorrEnergyCut(0.3);
299  }
300 
301 
302  // Create proper jet container
303  auto jetcont = task->AddJetContainer(jettype, AliJetContainer::antikt_algorithm, AliJetContainer::E_scheme, radius, acctype, tracks, clusters);
304  jetcont->SetName("datajets");
305  task->SetNameJetContainer("datajets");
306  std::cout << "Adding jet container with underlying array:" << jetcont->GetArrayName() << std::endl;
307 
308  // Link input and output container
309  outfilename << mgr->GetCommonFileName() << ":JetSpectrum_" << tag.str().data();
310  mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
311  mgr->ConnectOutput(task, 1, mgr->CreateContainer(Form("JetSpectrum_%s", tag.str().data()), TList::Class(), AliAnalysisManager::kOutputContainer, outfilename.str().data()));
312 
313  return task;
314 }
EMCal fiducial acceptance (each eta, phi edge narrowed by jet R)
Class creating a variable binning, used in the histogram manager.
Class creating a linear binning, used in the histogram manager.
AliJetContainer * GetJetContainer(Int_t i=0) const
Double_t Eta() const
Definition: AliEmcalJet.h:121
Double_t Phi() const
Definition: AliEmcalJet.h:117
Container with name, TClonesArray and cuts for particles.
void SetUseAliAnaUtils(Bool_t b, Bool_t bRejPilup=kTRUE)
Bool_t fRequireSubsetMB
Require for triggers to be a subset of Min. Bias (for efficiency studies)
static AliAnalysisTaskEmcalJetEnergySpectrum * AddTaskJetEnergySpectrum(Bool_t isMC, AliJetContainer::EJetType_t jettype, double radius, const char *trigger, const char *suffix="")
std::vector< TriggerCluster_t > GetTriggerClusterIndices(const TString &triggerstring) const
Double_t E() const
Definition: AliEmcalJet.h:119
Interface for binnings used by the histogram handler.
Definition: TBinning.h:21
static std::vector< PWG::EMCAL::Triggerinfo > DecodeTriggerString(const std::string &triggerstring)
Decoding trigger string.
void FillTHnSparse(const char *name, const double *x, double weight=1., Option_t *opt="")
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.
unsigned int UInt_t
Definition: External.C:33
THashList * GetListOfHistograms() const
Get the list of histograms.
Definition: THistManager.h:671
ULong_t fMinBiasTrigger
Min bias trigger for trigger subset (for efficiency studies)
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
Create a new TH1 within the container.
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
Bool_t fUseTriggerSelectionForData
Use trigger selection on data (require trigger patch in addition to trigger selection string) ...
Bool_t isMC
AliEmcalList * fOutput
!output list
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 test(int runnumber=195345)
void UserCreateOutputObjects()
Main initialization function on the worker.
bool Bool_t
Definition: External.C:53
Double_t NEF() const
Definition: AliEmcalJet.h:148
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)
static TString ClusterContainerNameFactory(Bool_t isAOD)
Get name of the default cluster container.
static TString TrackContainerNameFactory(Bool_t isAOD)
Get name of the default track container.