AliPhysics  c7b8e89 (c7b8e89)
AliAnalysisTaskEmcalJetHPerformance.cxx
Go to the documentation of this file.
1 
7 
8 #include <map>
9 #include <vector>
10 #include <iostream>
11 
12 #include <TObject.h>
13 #include <TCollection.h>
14 #include <TAxis.h>
15 #include <THnSparse.h>
16 
17 #include <AliAnalysisManager.h>
18 #include <AliLog.h>
19 
20 #include "yaml-cpp/yaml.h"
21 #include "THistManager.h"
22 #include "AliJetContainer.h"
23 #include "AliEmcalJet.h"
25 
27 
28 namespace PWGJE {
29 namespace EMCALJetTasks {
30 
35  AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetHPerformance", kFALSE),
36  fYAMLConfig(),
37  fConfigurationPath(""),
38  fConfigurationInitialized(false),
39  fHistManager(),
40  fEmbeddingQA(),
41  fCreateResponseMatrix(false),
42  fResponseMatrixFillMap(),
43  fResponseFromThreeJetCollections(true),
44  fMinFractionShared(0.),
45  fLeadingHadronBiasType(AliAnalysisTaskEmcalJetHUtils::kCharged)
46 {
47 }
48 
53  AliAnalysisTaskEmcalJet(name, kTRUE),
54  fYAMLConfig(),
57  fHistManager(name),
58  fCreateResponseMatrix(false),
63 {
64  // Ensure that additional general histograms are created
66 }
67 
72  fYAMLConfig(other.fYAMLConfig),
75  fHistManager(other.fHistManager.GetName()),
81 {
82  TIter next(other.fHistManager.GetListOfHistograms());
83  TObject* obj = 0;
84  while ((obj = next())) {
85  fHistManager.SetObject(obj, obj->GetName());
86  }
87 }
88 
94 {
95  swap(*this, other);
96  return *this;
97 }
98 
103 {
104  // Same ordering as in the constructor (for consistency)
105  std::string baseName = "enable";
106  fYAMLConfig.GetProperty({baseName, "responseMatrix"}, fCreateResponseMatrix, false);
107 
108  // Response matrix properties
109  baseName = "responseMatrix";
110  fYAMLConfig.GetProperty({baseName, "useThreeJetCollections"}, fResponseFromThreeJetCollections, false);
111  fYAMLConfig.GetProperty({baseName, "minFractionSharedPt"}, fMinFractionShared, false);
112  std::string hadronBiasStr = "";
113  baseName = "jets";
114  bool res = fYAMLConfig.GetProperty({baseName, "leadingHadronBiasType"}, hadronBiasStr, false);
115  // Only attempt to set the property if it is retrieved successfully
116  if (res) {
118  }
119 
120  // Base class options
121  // Recycle unused embedded events
122  fYAMLConfig.GetProperty("recycleUnusedEmbeddedEventsMode", fRecycleUnusedEmbeddedEventsMode, false);
123  // Centrality
124  std::pair<double, double> centRange;
125  res = fYAMLConfig.GetProperty("centralityRange", centRange, false);
126  if (res) {
127  SetCentRange(centRange.first, centRange.second);
128  }
129 }
130 
135 {
136  std::string baseName = "jets";
137  std::vector <std::string> jetNames = {"hybridLevelJets", "detLevelJets", "partLevelJets"};
138  for (const auto & jetName : jetNames) {
139  // Retrieve the node just to see if it is exists. If so, then we can proceed
140  YAML::Node node;
141  bool res = fYAMLConfig.GetProperty(std::vector<std::string>({baseName, jetName}), node, false);
142  if (res) {
143  // Retrieve jet container properties
144  std::string collectionName = "", acceptance = "";
145  double R = -1;
146  fYAMLConfig.GetProperty({baseName, jetName, "collection"}, collectionName, true);
147  fYAMLConfig.GetProperty({baseName, jetName, "acceptance"}, acceptance, true);
148  fYAMLConfig.GetProperty({baseName, jetName, "R"}, R, true);
149 
150  // Create jet container and set the name
151  AliDebugStream(1) << "Creating jet from jet collection name " << collectionName << " with acceptance " << acceptance << " and R=" << R << "\n";
152  AliJetContainer * jetCont = AddJetContainer(collectionName.c_str(), acceptance.c_str(), R);
153  jetCont->SetName(jetName.c_str());
154 
155  // Leading hadron type
156  int leadingHadronType = -1;
157  bool res = fYAMLConfig.GetProperty({baseName, jetName, "leadingHadronType"}, leadingHadronType, false);
158  if (res) {
159  AliDebugStream(1) << "Setting leading hadron type of " << leadingHadronType << " for jet cont " << jetName << "\n";
160  jetCont->SetLeadingHadronType(leadingHadronType);
161  }
162  }
163  else {
164  AliInfoStream() << "Unable to find definition of jet container corresponding to \"" << jetName << "\"\n";
165  }
166  }
167 }
168 
173 {
175 
176  // Setup and initialize the YAML config
177  if (fConfigurationPath != "") {
178  AliInfoStream() << "Adding YAML configuration found at \"" << fConfigurationPath << "\"\n";
179  int configPosition = fYAMLConfig.AddConfiguration(fConfigurationPath, "yamlConfig");
180  if (configPosition < 0) {
181  // Return immediately
183  }
184  }
185  else {
186  AliInfoStream() << "No YAML configuration fileanme passed.\n";
187  }
188 
189  // Always initialize for streaming purposes
191 
192  // Setup task based on the properties defined in the YAML config
193  AliDebugStream(2) << "Configuring task from the YAML configuration.\n";
196  AliDebugStream(2) << "Finished configuring via the YAML configuration.\n";
197 
198  // Print the results of the initialization
199  // Print outside of the ALICE Log system to ensure that it is always available!
200  std::cout << *this;
201 
204 }
205 
210 {
211  // Check that the task was initialized
213  AliFatal("Task was not initialized. Please ensure that Initialize() was called!");
214  }
215 
216  // Reinitialize the YAML configuration
218 
219  // Call the base class
221 
222  // Create histograms
223  if (fCreateResponseMatrix) {
225  }
226 
227  // Store hist manager output in the output list
228  TIter next(fHistManager.GetListOfHistograms());
229  TObject* obj = 0;
230  while ((obj = next())) {
231  fOutput->Add(obj);
232  }
233 
234  // Initialize
236  if (embeddingHelper) {
237  bool res = fEmbeddingQA.Initialize();
238  if (res) {
240  }
241  }
242 
243  // Post the data when finished
244  PostData(1, fOutput);
245 }
246 
251 {
252  // Main response matrix THnSparse
253  std::string name = "fHistResponseMatrix";
254  std::string title = "fHistResponseMatrix";
255  //histname = TString::Format("%s/BackgroundHistograms/hScaleFactorEMCal", jets->GetArrayName().Data());
256 
257  // Retrieve binning from the YAML configuration
258  std::vector<TAxis *> binning;
259  // This structure is designed to preserve the order of the axis in the YAML by using a YAML sequence (decoded into
260  // a vector), while defining a pair of the axis name and axis limts. Using this structure avoids the need to create
261  // a new object and conversion to retrieve the data
262  std::vector<std::pair<std::string, std::vector<double>>> sparseAxes;
263  std::string baseName = "responseMatrix";
264  fYAMLConfig.GetProperty({baseName, "axes"}, sparseAxes, true);
265  for (auto axis : sparseAxes) {
266  auto axisLimits = axis.second;
267  AliDebugStream(3) << "Creating axis " << axis.first << " with nBins " << axisLimits.at(0) << ", min: " << axisLimits.at(1) << ", max: " << axisLimits.at(2) << "\n";
268  binning.emplace_back(new TAxis(axisLimits.at(0), axisLimits.at(1), axisLimits.at(2)));
269  }
270 
271  // "s" ensures that Sumw2() is called
272  // The explicit const_cast is required
273  THnSparse * hist = fHistManager.CreateTHnSparse(name.c_str(), title.c_str(), binning.size(), const_cast<const TAxis **>(binning.data()), "s");
274  // Set the axis titles
275  int axisNumber = 0;
276  for (auto axis = sparseAxes.begin(); axis != sparseAxes.end(); axis++) {
277  AliDebugStream(5) << "ResponseMatrix: Add axis " << axis->first << " to sparse\n";
278  hist->GetAxis(axisNumber)->SetTitle(axis->first.c_str());
279  axisNumber++;
280  }
281 
282  // Define mapping from value name to value (including jet information)
283  for (unsigned int iJet = 1; iJet < 3; iJet++)
284  {
285  // pt
286  // ex: p_{T,1}
287  fResponseMatrixFillMap.insert(std::make_pair(std::string("p_{T,") + std::to_string(iJet) + "}", std::make_pair(iJet, &ResponseMatrixFillWrapper::fPt)));
288  // Area
289  // ex: A_{jet,1}
290  fResponseMatrixFillMap.insert(std::make_pair(std::string("A_{jet,") + std::to_string(iJet) + "}", std::make_pair(iJet, &ResponseMatrixFillWrapper::fArea)));
291  // EP angle
292  // ex: #theta_{jet,1}^{EP}
293  fResponseMatrixFillMap.insert(std::make_pair(std::string("#theta_{jet,") + std::to_string(iJet) + "}^{EP}", std::make_pair(iJet, &ResponseMatrixFillWrapper::fRelativeEPAngle)));
294  // Leading hadron
295  // ex: p_{T,particle,1}^{leading} (GeV/c)
296  fResponseMatrixFillMap.insert(std::make_pair(std::string("p_{T,particle,") + std::to_string(iJet) + "}^{leading} (GeV/c)", std::make_pair(iJet, &ResponseMatrixFillWrapper::fLeadingHadronPt)));
297  }
298  // Distance from one jet to another
299  fResponseMatrixFillMap.insert(std::make_pair("distance", std::make_pair(1, &ResponseMatrixFillWrapper::fDistance)) );
300  // Centrality
301  fResponseMatrixFillMap.insert(std::make_pair("centrality", std::make_pair(1, &ResponseMatrixFillWrapper::fCentrality)) );
302 
303  // Shared momentum fraction
304  name = "fHistFractionSharedPt";
305  title = "Fraction of p_{T} shared between matched jets";
306  // Need to include the bin from 1-1.01 to ensure that jets which shared all of their momentum
307  // due not end up in the overflow bin!
308  fHistManager.CreateTH1(name.c_str(), title.c_str(), 101, 0, 1.01);
309 }
310 
315 {
316  // Only fill the embedding qa plots if:
317  // - We are using the embedding helper
318  // - The class has been initialized
319  if (fEmbeddingQA.IsInitialized()) {
321  }
322 
323  if (fCreateResponseMatrix) {
325  }
326 
327  return kTRUE;
328 }
329 
334 {
335  AliJetContainer * jetsHybrid = GetJetContainer("hybridLevelJets");
336  AliJetContainer * jetsDetLevel = GetJetContainer("detLevelJets");
337  AliJetContainer * jetsPartLevel = GetJetContainer("partLevelJets");
338  if (!jetsHybrid) {
339  AliErrorStream() << "Could not retrieve hybrid jet collection.\n";
340  return;
341  }
342  if (!jetsDetLevel) {
343  AliErrorStream() << "Could not retrieve det level jet collection.\n";
344  return;
345  }
346  if (fResponseFromThreeJetCollections && !jetsPartLevel) {
347  AliErrorStream() << "Could not retrieve part level jet collection.\n";
348  return;
349  }
350 
351  // Handle matching of jets.
352  for (auto jet1 : jetsHybrid->accepted())
353  {
354  // Get jet the det level jet from the hybrid jet
355  AliEmcalJet * jet2 = jet1->ClosestJet();
356  if(!jet2) continue;
357 
358  AliDebugStream(4) << "jet2: " << jet2->toString() << "\n";
359 
360  // Check shared fraction
361  double sharedFraction = jetsHybrid->GetFractionSharedPt(jet1);
362  fHistManager.FillTH1("fHistFractionSharedPt", sharedFraction);
363  if (sharedFraction < fMinFractionShared) {
364  AliDebugStream(4) << "Rejecting jet due to momentum fraction of " << sharedFraction << ", smaller than the minimum.\n";
365  continue;
366  }
367  else {
368  AliDebugStream(4) << "Jet passed momentum fraction cut with value of " << sharedFraction << "\n";
369  }
370 
371  // Apply additional selection to jet 2
372  // TODO: Should we apply acceptance criteria to jet 2 here?
373 
374  // Get MC level jet
375  AliEmcalJet * jetToPass = 0;
377  AliEmcalJet * jet3 = jet2->ClosestJet();
378 
379  // Accept jet 3
380  UInt_t rejectionReason = 0;
381  if (!jetsPartLevel->AcceptJet(jet3, rejectionReason)) {
382  // TODO: Store rejection reasons
383  //fHistRejectionReason2->Fill(jets2->GetRejectionReasonBitPosition(rejectionReason), jet2->Pt());
384  continue;
385  }
386 
387  AliDebugStream(4) << "jet3: " << jet3->toString() << "\n";
388 
389  // Use for the response
390  AliDebugStream(4) << "Using part level jet for response\n";
391  jetToPass = jet3;
392  }
393  else {
394  // Use for the response
395  AliDebugStream(4) << "Using det level jet for response\n";
396  jetToPass = jet2;
397  }
398 
399  // Fill response
400  FillResponseMatrix(jet1, jetToPass);
401  }
402 
403 }
404 
409 {
410  if (!jet1 || !jet2) {
411  AliErrorStream() << "Null jet passed to fill response matrix";
412  }
413 
414  AliDebugStream(3) << "About to create ResponseMatrixFillWrappers\n";
415  AliDebugStream(4) << "jet1: " << jet1->toString() << "\n";
416  AliDebugStream(4) << "jet2: " << jet2->toString() << "\n";
417  // Create map from jetNumber to jet and initialize the objects
418  std::map<unsigned int, ResponseMatrixFillWrapper> jetNumberToJet = {
419  std::make_pair(1, CreateResponseMatrixFillWrapper(jet1)),
420  std::make_pair(2, CreateResponseMatrixFillWrapper(jet2))
421  };
422 
423  // Fill histograms
424  std::string histName = "fHistResponseMatrix";
425  std::vector<double> values;
426  THnSparse * response = static_cast<THnSparse*>(fHistManager.FindObject(histName.c_str()));
427  AliDebugStream(3) << "About to fill response matrix values\n";
428  AliDebugStream(4) << "jet1: " << jet1->toString() << "\n";
429  AliDebugStream(4) << "jet2: " << jet2->toString() << "\n";
430  for (unsigned int i = 0; i < response->GetNdimensions(); i++) {
431  std::string title = response->GetAxis(i)->GetTitle();
432 
433  // Retrieve pair of jet and pointer to extract the fill value
434  auto jetPair = fResponseMatrixFillMap.find(title);
435  if (jetPair != fResponseMatrixFillMap.end()) {
436  auto wrapper = jetNumberToJet.at(jetPair->second.first);
437  auto member = jetPair->second.second;
438  AliDebugStream(4) << "Filling value " << wrapper.*member << " into axis " << title << "\n";
439  values.emplace_back(wrapper.*member);
440  }
441  else {
442  AliWarningStream() << "Unable to fill dimension " << title << "!\n";
443  }
444  }
445 
446  fHistManager.FillTHnSparse(histName.c_str(), values.data());
447 }
448 
453 {
455  if (!jet) {
456  AliErrorStream() << "Must pass valid jet to create object.\n";
457  return wrapper;
458  }
459  wrapper.fPt = jet->Pt();
460  wrapper.fArea = jet->Area();
461  wrapper.fPhi = jet->Phi();
462  wrapper.fDistance = jet->ClosestJetDistance();
463  wrapper.fCentrality = fCent;
466 
467  return wrapper;
468 }
469 
476 {
477  // Get the pointer to the existing analysis manager via the static access method.
478  //==============================================================================
479  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
480  if (!mgr)
481  {
482  AliErrorClass("No analysis manager to connect to.");
483  return nullptr;
484  }
485 
486  // Setup task name
487  std::string taskName = "AliAnalysisTaskEmcalJetHPerformance";
488  std::string suffixName(suffix);
489  if (suffixName != "") {
490  taskName += "_";
491  taskName += suffixName;
492  }
493 
494  // Create task and configure as desired
496  // Configure via the YAML configuration
497  mgr->AddTask(task);
498 
499  // Create containers for input/output
500  mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer() );
501  AliAnalysisDataContainer * outputContainer = mgr->CreateContainer(task->GetName(),
502  TList::Class(),
503  AliAnalysisManager::kOutputContainer,
504  Form("%s", AliAnalysisManager::GetCommonFileName()));
505  mgr->ConnectOutput(task, 1, outputContainer);
506 
507  return task;
508 }
509 
516 {
517  std::stringstream tempSS;
518  tempSS << std::boolalpha;
519  tempSS << "Jet collections:\n";
520  TIter next(&fJetCollArray);
521  AliJetContainer * jetCont;
522  while ((jetCont = static_cast<AliJetContainer *>(next()))) {
523  tempSS << "\t" << jetCont->GetName() << ": " << jetCont->GetArrayName() << "\n";
524  }
525  tempSS << "Response matrix:\n";
526  tempSS << "\tEnabled: " << fCreateResponseMatrix << "\n";
527  tempSS << "\tConstruct response from 3 jet collections: " << fResponseFromThreeJetCollections << "\n";
528  tempSS << "\tMin fraction shared pt: " << fMinFractionShared << "\n";
529  tempSS << "\tJet leading hadron bias type: " << fLeadingHadronBiasType << "\n";
530  tempSS << "\tResponse matrix fill map: \n";
531  for (auto el : fResponseMatrixFillMap) {
532  tempSS << "\t\tProperty " << el.first << " applied to jet " << el.second.first << "\n";
533  }
534 
535  return tempSS.str();
536 }
537 
544 std::ostream & AliAnalysisTaskEmcalJetHPerformance::Print(std::ostream & in) const {
545  in << toString();
546  return in;
547 }
548 
556 {
557  Printf("%s", toString().c_str());
558 }
559 
560 } /* namespace EMCALJetTasks */
561 } /* namespace PWGJE */
562 
570 std::ostream & operator<<(std::ostream & in, const PWGJE::EMCALJetTasks::AliAnalysisTaskEmcalJetHPerformance & myTask)
571 {
572  std::ostream & result = myTask.Print(in);
573  return result;
574 }
575 
580 {
581  using std::swap;
582 
583  // Same ordering as in the constructors (for consistency)
584  swap(first.fYAMLConfig, second.fYAMLConfig);
587  //swap(first.fHistManager, second.fHistManager); // Skip here, as there is no copy constructor for THistManager
593 }
594 
void SetCentRange(Double_t min, Double_t max)
Double_t Area() const
Definition: AliEmcalJet.h:130
bool fResponseFromThreeJetCollections
If true, the det level jets in collection 2 are only an intermediate step. They are used to get part ...
const char * title
Definition: MakeQAPdf.C:27
AliEmcalJet * ClosestJet() const
Definition: AliEmcalJet.h:327
AliJetContainer * GetJetContainer(Int_t i=0) const
std::map< std::string, std::pair< int, double ResponseMatrixFillWrapper::* > > fResponseMatrixFillMap
! Map from axis title to pair of (jet number, function to retrieve fill value)
Double_t ClosestJetDistance() const
Definition: AliEmcalJet.h:328
Bool_t fRecycleUnusedEmbeddedEventsMode
Allows the recycling of embedded events which fail internal event selection. See the embedding helper...
void SetLeadingHadronType(Int_t t)
Double_t Phi() const
Definition: AliEmcalJet.h:117
double fCentrality
! Centrality of the given jet (an event level property, but useful to here available here) ...
Double_t fEPV0
!event plane V0
Declaration of class AliAnalysisTaskEmcalEmbeddingHelper.
std::string fConfigurationPath
Path to the YAML configuration file.
virtual Bool_t AcceptJet(Int_t i, UInt_t &rejectionReason) const
AliJetContainer * AddJetContainer(const char *n, TString defaultCutType, Float_t jetRadius=0.4)
AliAnalysisTaskEmcalJetHPerformance & operator=(AliAnalysisTaskEmcalJetHPerformance other)
bool GetProperty(std::vector< std::string > propertyPath, const std::string &propertyName, T &property, const bool requiredProperty) const
AliEmcalEmbeddingQA fEmbeddingQA
! Embedding QA hists (will only be added if embedding)
void SetObject(TObject *const o, const char *group="/")
Set a new group into the container into the parent group.
ResponseMatrixFillWrapper CreateResponseMatrixFillWrapper(AliEmcalJet *jet) const
PWG::Tools::AliYAMLConfiguration fYAMLConfig
YAML configuration file.
bool fCreateResponseMatrix
If true, create a response matrix with the available jet collections.
bool AddQAPlotsToList(TList *list)
AliAnalysisTaskEmcalJetHUtils::ELeadingHadronBiasType_t fLeadingHadronBiasType
Leading hadron in jet bias type (either charged, neutral, or both)
void FillTHnSparse(const char *name, const double *x, double weight=1., Option_t *opt="")
bool IsInitialized() const
TObject * FindObject(const char *name) const
Find an object inside the container.
bool fConfigurationInitialized
True if the task configuration has been successfully initialized.
unsigned int UInt_t
Definition: External.C:33
THashList * GetListOfHistograms() const
Get the list of histograms.
Definition: THistManager.h:671
Implementation of task to embed external events.
TString toString() const
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
Create a new TH1 within the container.
Double_t fCent
!event centrality
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
Fill a 1D histogram within the container.
static AliAnalysisTaskEmcalJetHPerformance * AddTaskEmcalJetHPerformance(const char *suffix="")
static double RelativeEPAngle(double jetAngle, double epAngle)
TObjArray fJetCollArray
jet collection array
static double GetLeadingHadronPt(AliEmcalJet *jet, ELeadingHadronBiasType_t leadingHadronType)
Double_t Pt() const
Definition: AliEmcalJet.h:109
AliEmcalList * fOutput
!output list
void SetMakeGeneralHistograms(Bool_t g)
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:51
int AddConfiguration(std::string configurationFilename, std::string configurationName="")
Double_t GetFractionSharedPt(const AliEmcalJet *jet, AliParticleContainer *cont2=0x0) const
const char Option_t
Definition: External.C:48
void UserCreateOutputObjects()
Main initialization function on the worker.
const AliJetIterableContainer accepted() const
bool Bool_t
Definition: External.C:53
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.
void swap(PWGJE::EMCALJetTasks::AliAnalysisTaskEmcalJetHPerformance &first, PWGJE::EMCALJetTasks::AliAnalysisTaskEmcalJetHPerformance &second)
double fMinFractionShared
Minimum fraction of shared jet pt required for matching a hybrid jet to detector level.
Container for jet within the EMCAL jet framework.
static const std::map< std::string, ELeadingHadronBiasType_t > fgkLeadingHadronBiasMap
! Map from name to leading hadron bias used with the YAML config
friend std::ostream & operator<<(std::ostream &in, const AliAnalysisTaskEmcalJetHPerformance &myTask)
static const AliAnalysisTaskEmcalEmbeddingHelper * GetInstance()