AliPhysics  master (3d17d9d)
AliAnalysisTaskEmcalJetHPerformance.cxx
Go to the documentation of this file.
1 
6 
7 #include <cmath>
8 #include <map>
9 #include <vector>
10 #include <iostream>
11 #include <bitset>
12 
13 #include <TObject.h>
14 #include <TCollection.h>
15 #include <TAxis.h>
16 #include <THnSparse.h>
17 
18 #include <AliAnalysisManager.h>
19 #include <AliInputEventHandler.h>
20 #include <AliLog.h>
21 
22 #include "yaml-cpp/yaml.h"
24 #include "AliTrackContainer.h"
25 #include "AliTLorentzVector.h"
26 #include "AliClusterContainer.h"
27 #include "AliEmcalContainerUtils.h"
28 #include "AliJetContainer.h"
29 #include "AliEmcalJet.h"
30 #include "AliAnalysisTaskFlowVectorCorrections.h"
31 #include "AliQnCorrectionsManager.h"
32 #include "AliQnCorrectionsQnVector.h"
33 
35 
37 
38 namespace PWGJE {
39 namespace EMCALJetTasks {
40 
45  AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetHPerformance", kFALSE),
46  fYAMLConfig(),
47  fConfigurationInitialized(false),
48  fHistManager(),
49  fEmbeddingQA(),
50  fCreateQAHists(false),
51  fCreateCellQAHists(false),
52  fPerformJetMatching(false),
53  fCreateResponseMatrix(false),
54  fEmbeddedCellsName("emcalCells"),
55  fPreviousEventTrigger(0),
56  fPreviousEmbeddedEventSelected(false),
57  fEfficiencyPeriodIdentifier(AliAnalysisTaskEmcalJetHUtils::kDisableEff),
58  fFlowQnVectorManager(nullptr),
59  fResponseMatrixFillMap(),
60  fResponseFromThreeJetCollections(true),
61  fMinFractionShared(0.),
62  fLeadingHadronBiasType(AliAnalysisTaskEmcalJetHUtils::kCharged)
63 {
64 }
65 
70  AliAnalysisTaskEmcalJet(name, kTRUE),
71  fYAMLConfig(),
73  fHistManager(name),
74  fEmbeddingQA(),
75  fCreateQAHists(false),
76  fCreateCellQAHists(false),
77  fPerformJetMatching(false),
78  fCreateResponseMatrix(false),
79  fEmbeddedCellsName("emcalCells"),
88 {
89  // Ensure that additional general histograms are created
91 }
92 
97  fYAMLConfig(other.fYAMLConfig),
99  fHistManager(other.fHistManager.GetName()),
100  //fEmbeddingQA(other.fEmbeddingQA), // Cannot use because the THistManager (which is in the class) copy constructor is private.
114 {
115  TIter next(other.fHistManager.GetListOfHistograms());
116  TObject* obj = 0;
117  while ((obj = next())) {
118  fHistManager.SetObject(obj, obj->GetName());
119  }
120 }
121 
127 {
128  swap(*this, other);
129  return *this;
130 }
131 
141 {
143  trackPt, trackEta, fCentBin, fEfficiencyPeriodIdentifier,
144  "PWGJE::EMCALJetTasks::AliAnalysisTaskEmcalJetHPerformance");
145 }
146 
151 {
152  // Base class options
153  // Task physics (trigger) selection.
154  std::vector<std::string> physicsSelection;
155  bool res = fYAMLConfig.GetProperty(std::vector<std::string>({"eventCuts", "physicsSelection"}), physicsSelection, false);
156  if (res) {
157  fOfflineTriggerMask = AliEmcalContainerUtils::DeterminePhysicsSelectionFromYAML(physicsSelection);
158  }
159 
160  // Same ordering as in the constructor (for consistency)
161  std::string baseName = "enable";
162  // These are disabled by default.
163  fYAMLConfig.GetProperty({baseName, "QAHists"}, fCreateQAHists, false);
164  fYAMLConfig.GetProperty({baseName, "CellQAHists"}, fCreateCellQAHists, false);
165  fYAMLConfig.GetProperty({baseName, "jetMatching"}, fPerformJetMatching, false);
166  fYAMLConfig.GetProperty({baseName, "responseMatrix"}, fCreateResponseMatrix, false);
167 
168  // Event cuts
169  baseName = "eventCuts";
170  // If event cuts are enabled (which they exceptionally are by default), then we want to configure them here.
171  // If the event cuts are explicitly disabled, then we invert that value to enable the AliAnylsisTaskEmcal
172  // builtin event selection.
173  bool tempBool;
174  fYAMLConfig.GetProperty({baseName, "enabled"}, tempBool, false);
175  fUseBuiltinEventSelection = !tempBool;
176  if (fUseBuiltinEventSelection == false) {
177  // Need to include the namespace so that AliDebug will work properly...
178  std::string taskName = "PWGJE::EMCALJetTasks::";
179  taskName += GetName();
180  AliAnalysisTaskEmcalJetHUtils::ConfigureEventCuts(fAliEventCuts, fYAMLConfig, fOfflineTriggerMask, baseName, taskName);
181  }
182 
183  // General task options
184  baseName = "general";
185  fYAMLConfig.GetProperty({baseName, "nCentBins"}, fNcentBins, false);
186 
187  // QA options
188  baseName = "QA";
189  fYAMLConfig.GetProperty({baseName, "cellsName"}, fCaloCellsName, false);
190  // Defaults to "emcalCells" if not set.
191  fYAMLConfig.GetProperty({baseName, "embeddedCellsName"}, fEmbeddedCellsName, false);
192  // Efficiency
193  std::string tempStr = "";
194  baseName = "efficiency";
195  res = fYAMLConfig.GetProperty({baseName, "periodIdentifier"}, tempStr, false);
196  if (res) {
198  }
199 
200  // Jet matching
201  baseName = "jetMatching";
202  fYAMLConfig.GetProperty({baseName, "maxMatchingDistance"}, fMaxJetMatchingDistance, false);
203 
204  // Response matrix properties
205  baseName = "responseMatrix";
206  fYAMLConfig.GetProperty({baseName, "useThreeJetCollections"}, fResponseFromThreeJetCollections, false);
207  fYAMLConfig.GetProperty({baseName, "minFractionSharedPt"}, fMinFractionShared, false);
208  std::string hadronBiasStr = "";
209  baseName = "jets";
210  res = fYAMLConfig.GetProperty({baseName, "leadingHadronBiasType"}, hadronBiasStr, false);
211  // Only attempt to set the property if it is retrieved successfully
212  if (res) {
214  }
215 }
216 
221 {
222  std::string baseName = "jets";
223  std::vector <std::string> jetNames = {"hybridLevelJets", "detLevelJets", "partLevelJets", "analysisJets"};
224  for (const auto & jetName : jetNames) {
225  // Retrieve the node just to see if it is exists. If so, then we can proceed
226  YAML::Node node;
227  bool res = fYAMLConfig.GetProperty(std::vector<std::string>({baseName, jetName}), node, false);
228  if (res) {
229  // Retrieve jet container properties
230  std::string collectionName = "";
231  double R = -1;
232  fYAMLConfig.GetProperty({baseName, jetName, "collection"}, collectionName, true);
233  fYAMLConfig.GetProperty({baseName, jetName, "R"}, R, true);
234 
235  // Determine the jet acceptance.
236  std::vector<std::string> jetAcceptances;
237  fYAMLConfig.GetProperty({baseName, jetName, "acceptance"}, jetAcceptances, true);
239 
240  // Create jet container and set the name
241  AliDebugStream(1) << "Creating jet from jet collection name " << collectionName << " with acceptance "
242  << jetAcceptance << " and R=" << R << "\n";
243  AliJetContainer* jetCont = AddJetContainer(collectionName.c_str(), jetAcceptance, R);
244  jetCont->SetName(jetName.c_str());
245 
246  // Jet area cut percentage
247  double jetAreaCut = -1;
248  bool res = fYAMLConfig.GetProperty({ baseName, jetName, "areaCutPercentage" }, jetAreaCut, false);
249  if (res) {
250  AliDebugStream(1) << "Setting jet area cut percentage of " << jetAreaCut << " for jet cont " << jetName
251  << "\n";
252  jetCont->SetPercAreaCut(jetAreaCut);
253  }
254 
255  // Rho name
256  std::string tempStr = "";
257  res = fYAMLConfig.GetProperty({ baseName, jetName, "rhoName" }, tempStr, false);
258  if (res) {
259  AliDebugStream(1) << "Setting rho name of " << tempStr << " for jet cont " << jetName << "\n";
260  jetCont->SetRhoName(tempStr.c_str());
261  }
262 
263  // Leading hadron type
264  int leadingHadronType = -1;
265  res = fYAMLConfig.GetProperty({ baseName, jetName, "leadingHadronType" }, leadingHadronType, false);
266  if (res) {
267  AliDebugStream(1) << "Setting leading hadron type of " << leadingHadronType << " for jet cont "
268  << jetName << "\n";
269  jetCont->SetLeadingHadronType(leadingHadronType);
270  }
271 
272  // Leading hadron cut
273  double leadingHadronBias = 0.;
274  res = fYAMLConfig.GetProperty({ baseName, jetName, "leadingHadronBias" }, leadingHadronBias, false);
275  if (res) {
276  AliDebugStream(1) << "Setting leading hadron bias of " << leadingHadronBias << " for jet cont "
277  << jetName << "\n";
278  jetCont->SetMinTrackPt(leadingHadronBias);
279  }
280 
281  // Max track pt
282  double maxTrackPt = 0.;
283  res = fYAMLConfig.GetProperty({ baseName, jetName, "maxTrackPt" }, maxTrackPt, false);
284  if (res) {
285  AliDebugStream(1) << "Setting max track pt of " << maxTrackPt << " for jet cont "
286  << jetName << "\n";
287  jetCont->SetMinTrackPt(maxTrackPt);
288  }
289  }
290  else {
291  AliInfoStream() << "Unable to find definition of jet container corresponding to \"" << jetName << "\"\n";
292  }
293  }
294 }
295 
300 {
301  std::string baseName = "particles";
302  // Retrieve the node just to see if it is exists. If so, then we can proceed
303  YAML::Node node;
304  fYAMLConfig.GetProperty(baseName, node);
305  // Iterate over all of the particle and track containers
306  for (const auto & n : node) {
307  std::string containerName = n.first.as<std::string>();
308 
309  // Create the container.
310  std::string branchName = "";
311  fYAMLConfig.GetProperty({baseName, containerName, "branchName"}, branchName, true);
312  if (branchName == "usedefault") {
314  }
316  this->AdoptParticleContainer(partCont);
317 
318  // Configure the container
319  // Need to include the namespace so that AliDebug will work properly...
320  std::string taskName = "PWGJE::EMCALJetTasks::";
321  taskName += GetName();
322  std::vector<std::string> baseNameWithContainer = { baseName, containerName };
324  partCont, fYAMLConfig, taskName);
325 
326  // Track specific properties
327  AliTrackContainer * trackCont = dynamic_cast<AliTrackContainer *>(partCont);
328  if (trackCont) {
330  fYAMLConfig, taskName);
331  }
332 
333  AliDebugStream(2) << "Particle/track container: " << partCont->GetName()
334  << ", array class: " << partCont->GetClassName() << ", collection name: \""
335  << partCont->GetArrayName() << "\", min pt: " << partCont->GetMinPt() << "\n";
336  }
337 }
338 
343 {
344  std::string baseName = "clusters";
345  // Retrieve the node just to see if it is exists. If so, then we can proceed
346  YAML::Node node;
347  fYAMLConfig.GetProperty(baseName, node);
348  // Iterate over all of the cluster containers
349  for (const auto & n : node) {
350  std::string containerName = n.first.as<std::string>();
351 
352  // Create the container.
353  std::string branchName = "";
354  fYAMLConfig.GetProperty({baseName, containerName, "branchName"}, branchName, true);
355  if (branchName == "usedefault") {
357  }
358  AliClusterContainer * clusterCont = this->AddClusterContainer(branchName.c_str());
359 
360  // Configure the container
361  // Need to include the namespace so that AliDebug will work properly...
362  std::string taskName = "PWGJE::EMCALJetTasks::";
363  taskName += GetName();
364  std::vector<std::string> baseNameWithContainer = { baseName, containerName };
366  clusterCont, fYAMLConfig, taskName);
367 
368  // Cluster specific properties
369  AliAnalysisTaskEmcalJetHUtils::ConfigureClusterContainersFromYAMLConfig(baseNameWithContainer, clusterCont, fYAMLConfig, taskName);
370 
371  AliDebugStream(2) << "Cluster container: " << clusterCont->GetName()
372  << ", array class: " << clusterCont->GetClassName() << ", collection name: \""
373  << clusterCont->GetArrayName() << "\", min pt: " << clusterCont->GetMinPt() << "\n";
374  }
375 }
376 
381 {
383 
384  // Ensure that we have at least one configuration in the YAML config.
385  if (fYAMLConfig.DoesConfigurationExist(0) == false) {
386  // No configurations exist. Return immediately.
388  }
389 
390  // Always initialize for streaming purposes
392 
393  // Setup task based on the properties defined in the YAML config
394  AliDebugStream(2) << "Configuring task from the YAML configuration.\n";
399  AliDebugStream(2) << "Finished configuring via the YAML configuration.\n";
400 
401  // Print the results of the initialization
402  // Print outside of the ALICE Log system to ensure that it is always available!
403  std::cout << *this;
404 
407 }
408 
413 {
414  // First call the base class
416 
417  // Check that the task was initialized
419  AliFatal("Task was not initialized. Please ensure that Initialize() was called!");
420  }
421  // Reinitialize the YAML configuration
423 
424  // Create histograms
425  if (fCreateQAHists) {
426  SetupQAHists();
427  }
428  if (fPerformJetMatching) {
430  }
431  if (fCreateResponseMatrix) {
433  }
434 
435  // Store hist manager output in the output list
436  TIter next(fHistManager.GetListOfHistograms());
437  TObject* obj = 0;
438  while ((obj = next())) {
439  fOutput->Add(obj);
440  }
441 
442  // Initialize
444  if (embeddingHelper) {
445  bool res = fEmbeddingQA.Initialize();
446  if (res) {
448  }
449  }
450 
451  // Post the data when finished
452  PostData(1, fOutput);
453 }
454 
462 {
463  std::string name = prefix + "/fHistCellEnergy";
464  std::string title = name + ";\\mathit{E}_{\\mathrm{cell}} (\\mathrm{GeV});counts";
465  fHistManager.CreateTH1(name.c_str(), title.c_str(), 300, 0, 150);
466 
467  name = prefix + "/fHistCellTime";
468  title = name + ";t (s);counts";
469  fHistManager.CreateTH1(name.c_str(), title.c_str(), 1000, -10e-6, 10e-6);
470 
471  name = prefix + "/fHistNCells";
472  title = name + ";\\mathrm{N}_{\\mathrm{cells}};counts";
473  fHistManager.CreateTH1(name.c_str(), title.c_str(), 100, 0, 5000);
474 
475  name = prefix + "/fHistCellID";
476  title = name + ";\\mathrm{N}_{\\mathrm{cell}};counts";
477  fHistManager.CreateTH1(name.c_str(), title.c_str(), 20000, 0, 20000);
478 
479  // Histograms for embedding QA which use the cell timing to determine whether the
480  // embedded event has been double corrected.
481  if (prefix.find("embedding") != std::string::npos) {
482  name = prefix + "/fHistEmbeddedEventUsed";
483  title = name + ";Embedded event used";
484  fHistManager.CreateTH1(name.c_str(), title.c_str(), 2, 0, 2);
485 
486  name = prefix + "/fHistInternalEventSelection";
487  title = name + ";Embedded event used;Trigger bit";
488  fHistManager.CreateTH2(name.c_str(), title.c_str(), 2, 0, 2, 32, -0.5, 31.5);
489  }
490 }
491 
496 {
497  // Only create and fill the cells QA if we've setup the cells name.
498  if (fCaloCellsName != "") {
499  std::string prefix = "QA/";
500  prefix += fCaloCellsName.Data();
502  auto embeddingInstance = AliAnalysisTaskEmcalEmbeddingHelper::GetInstance();
503  if (embeddingInstance) {
504  prefix = "QA/embedding/";
505  prefix += fEmbeddedCellsName;
507  }
508  }
509 }
510 
515 {
516  // Cell level QA
517  if (fCreateCellQAHists) {
519  }
520 
521  // Tracks
522  AliTrackContainer * trackCont = nullptr;
523  TIter nextTrackColl(&fParticleCollArray);
524  while ((trackCont = static_cast<AliTrackContainer*>(nextTrackColl()))) {
525  std::string name = "QA/%s/fHistTrackPtEtaPhi";
526  std::string title = name + ";#it{p}_{T} (GeV);#eta;#phi";
527  fHistManager.CreateTH3(TString::Format(name.c_str(), trackCont->GetName()),
528  TString::Format(title.c_str(), trackCont-> GetName()), 50, 0, 25, 40, -1, 1, 72, 0,
529  TMath::TwoPi());
530  name = "QA/%s/fHistTrackPtEtaPhiEfficiencyCorrected";
531  title = name + ";#it{p}_{T} (GeV);#eta;#phi";
532  fHistManager.CreateTH3(TString::Format(name.c_str(), trackCont->GetName()),
533  TString::Format(title.c_str(), trackCont-> GetName()), 50, 0, 25, 40, -1, 1, 72, 0,
534  TMath::TwoPi());
535  }
536 
537  // Clusters
538  AliClusterContainer * clusterCont = nullptr;
539  TIter nextClusColl(&fClusterCollArray);
540  while ((clusterCont = static_cast<AliClusterContainer*>(nextClusColl()))) {
541  // Cluster time vs energy
542  std::string name = "QA/%s/fHistClusterEnergyVsTime";
543  std::string title = name + ";#it{E}_{cluster} (GeV);t_{cluster} (s);counts";
544  fHistManager.CreateTH2(TString::Format(name.c_str(), clusterCont->GetName()),
545  TString::Format(title.c_str(), clusterCont->GetName()), 1000, 0, 100, 300, -300e-9,
546  300e-9);
547  // Hadronically corrected energy (which is usually what we're using)
548  name = "QA/%s/fHistClusterHadCorrEnergy";
549  title = name + ";#it{E}_{cluster}^{had.corr.} (GeV);counts";
550  fHistManager.CreateTH1(TString::Format(name.c_str(), clusterCont->GetName()), TString::Format(title.c_str(), clusterCont->GetName()), 200, 0, 100);
551  // Cluster eta-phi
552  name = "QA/%s/fHistClusterEtaPhi";
553  title = name + ";#eta_{cluster};#phi_{cluster};counts";
554  fHistManager.CreateTH2(TString::Format(name.c_str(), clusterCont->GetName()), TString::Format(title.c_str(), clusterCont->GetName()), 28, -0.7, 0.7, 72, 0, TMath::TwoPi());
555  }
556 
557  // Jets
558  AliJetContainer * jetCont = nullptr;
559  TIter nextJetColl(&fJetCollArray);
560  while ((jetCont = static_cast<AliJetContainer*>(nextJetColl()))) {
561  // Jet pT
562  std::string name = "QA/%s/fHistJetPt";
563  std::string title = name + ";p_{T} (GeV)";
564  fHistManager.CreateTH1(TString::Format(name.c_str(), jetCont->GetName()), TString::Format(title.c_str(), jetCont->GetName()), 600, -50, 250);
565  }
566 
567  // Event plane resolution
568  // Only enable if the Qn vector corrections task is available.
569  auto qnVectorTask = dynamic_cast<AliAnalysisTaskFlowVectorCorrections*>(
570  AliAnalysisManager::GetAnalysisManager()->GetTask("FlowQnVectorCorrections"));
571  if (qnVectorTask) {
572  // Setup the Qn corrections manager.
573  fFlowQnVectorManager = qnVectorTask->GetAliQnCorrectionsManager();
574  // Then setup the resolution histograms.
575  // Resolution are calculated via TProfiles as a function of centrality.
576  std::map<std::string, std::string> resolutionNamesToTitles = {
577  std::make_pair("VZERO", R"(VZERO: $<\cos(R * (\Psi_{\mathrm{TPC_{Pos}}} - \Psi_{\mathrm{TPC_{Neg}}}))>$)"),
578  std::make_pair("TPC_Pos", R"(TPC $\eta$ > 0: $<\cos(R * (\Psi_{\mathrm{V0M}} - \Psi_{\mathrm{TPC_{Neg}}}))>$)"),
579  std::make_pair("TPC_Neg", R"(TPC $\eta$ < 0: $<\cos(R * (\Psi_{\mathrm{V0M}} - \Psi_{\mathrm{TPC_{Pos}}}))>$)")
580  };
581  for (const auto& nameAndTitle: resolutionNamesToTitles) {
582  for (unsigned int R = 2; R < 9; R++) {
583  std::string name = "QA/eventPlaneRes/%s/R%i";
584  std::string title = "%s event plane res;Centrality (%%);R_{%i}";
585  fHistManager.CreateTProfile(TString::Format(name.c_str(), nameAndTitle.first.c_str(), R),
586  TString::Format(title.c_str(), nameAndTitle.second.c_str(), R),
587  10, 0, 100);
588  }
589  }
590  }
591 }
592 
597 {
598  // Setup
599  const int nBinsPt = 40;
600  const int nBinsDPhi = 72;
601  const int nBinsDEta = 100;
602  const int nBinsDR = 50;
603  const int nBinsFraction = 101;
604 
605  const double minPt = -50.;
606  const double maxPt = 150.;
607  const double minDPhi = -0.5;
608  const double maxDPhi = 0.5;
609  const double minDEta = -0.5;
610  const double maxDEta = 0.5;
611  const double minDR = 0.;
612  const double maxDR = 0.5;
613  const double minFraction = -0.005;
614  const double maxFraction = 1.005;
615 
616  // Create the histograms
617  // "s" ensures that Sumw2() is called
618  std::string name = "";
619  std::string title = "";
620  std::vector<std::string> matchingTypes = {"hybridToDet", "detToPart"};
621  for (auto matchingType : matchingTypes)
622  {
623  for (unsigned int centBin = 0; centBin < fNcentBins; centBin++)
624  {
625  // Jet 1 pt vs delta eta vs delta phi
626  name = "jetMatching/%s/fh3PtJet1VsDeltaEtaDeltaPhi_%d";
627  title = "fh3PtJet1VsDeltaEtaDeltaPhi_%d;#it{p}_{T,jet1};#it{#Delta#eta};#it{#Delta#varphi}";
628  fHistManager.CreateTH3(TString::Format(name.c_str(), matchingType.c_str(), centBin),
629  TString::Format(title.c_str(), centBin), nBinsPt, minPt, maxPt, nBinsDEta,
630  minDEta, maxDEta, nBinsDPhi, minDPhi, maxDPhi, "s");
631 
632  // Jet pt 1 vs delta R
633  name = "jetMatching/%s/fh2PtJet1VsDeltaR_%d";
634  title = "fh2PtJet1VsDeltaR_%d;#it{p}_{T,jet1};#it{#Delta R}";
635  fHistManager.CreateTH2(TString::Format(name.c_str(), matchingType.c_str(), centBin),
636  TString::Format(title.c_str(), centBin), nBinsPt, minPt, maxPt, nBinsDR, minDR,
637  maxDR, "s");
638 
639  // Jet pt 2 vs shared momentum fraction
640  name = "jetMatching/%s/fh2PtJet2VsFraction_%d";
641  title = "fh2PtJet2VsFraction_%d;#it{p}_{T,jet2};#it{f}_{shared}";
642  fHistManager.CreateTH2(TString::Format(name.c_str(), matchingType.c_str(), centBin),
643  TString::Format(title.c_str(), centBin), nBinsPt, minPt, maxPt, nBinsFraction,
644  minFraction, maxFraction, "s");
645 
646  // Jet pt 1 vs lead pt before checking if the jet has a matched jet.
647  name = "jetMatching/%s/fh2PtJet1VsLeadPtAllSel_%d";
648  title = "fh2PtJet1VsLeadPtAllSel_%d;#it{p}_{T,jet1};#it{p}_{T,lead trk}";
649  fHistManager.CreateTH2(TString::Format(name.c_str(), matchingType.c_str(), centBin),
650  TString::Format(title.c_str(), centBin), nBinsPt, minPt, maxPt, 20, 0., 20.,
651  "s");
652 
653  // Jet pt 1 vs lead pt after checking if the jet has a matched jet.
654  name = "jetMatching/%s/fh2PtJet1VsLeadPtTagged_%d";
655  title = "fh2PtJet1VsLeadPtTagged_%d;#it{p}_{T,jet1};#it{p}_{T,lead trk}";
656  fHistManager.CreateTH2(TString::Format(name.c_str(), matchingType.c_str(), centBin),
657  TString::Format(title.c_str(), centBin), nBinsPt, minPt, maxPt, 20, 0., 20.,
658  "s");
659 
660  // Jet pt 1 vs jet pt 2.
661  name = "jetMatching/%s/fh2PtJet1VsPtJet2_%d";
662  title = "fh2PtJet1VsPtJet2_%d;#it{p}_{T,jet1};#it{p}_{T,jet2}";
663  fHistManager.CreateTH2(TString::Format(name.c_str(), matchingType.c_str(), centBin),
664  TString::Format(title.c_str(), centBin), nBinsPt, minPt, maxPt, nBinsPt, minPt,
665  maxPt, "s");
666 
667  // Jet pt2 residual (aka. the jet energy scale).
668  name = "jetMatching/%s/fh2PtJet2VsRelPt_%d";
669  title = "fh2PtJet2VsRelPt_%d;#it{p}_{T,jet2};(#it{p}_{T,jet1}-#it{p}_{T,jet2})/#it{p}_{T,jet2}";
670  fHistManager.CreateTH2(TString::Format(name.c_str(), matchingType.c_str(), centBin),
671  TString::Format(title.c_str(), centBin), nBinsPt, minPt, maxPt, 241, -2.41, 2.41,
672  "s");
673  }
674 
675  // Number of jets accepted.
676  name = "jetMatching/%s/fNAccJets";
677  title = "fNAccJets;N/ev";
678  fHistManager.CreateTH1(TString::Format(name.c_str(), matchingType.c_str()), title.c_str(), 11, -0.5, 9.5, "s");
679  }
680 
681 }
682 
687 {
688  // Main response matrix THnSparse
689  std::string name = "response/fHistResponseMatrix";
690  std::string title = name;
691 
692  // Retrieve binning from the YAML configuration
693  std::vector<TAxis *> binning;
694  // This structure is designed to preserve the order of the axis in the YAML by using a YAML sequence (decoded into
695  // a vector), while defining a pair of the axis name and axis limts. Using this structure avoids the need to create
696  // a new object and conversion to retrieve the data
697  std::vector<std::pair<std::string, std::vector<double>>> sparseAxes;
698  std::string baseName = "responseMatrix";
699  fYAMLConfig.GetProperty({baseName, "axes"}, sparseAxes, true);
700  for (auto axis : sparseAxes) {
701  auto axisLimits = axis.second;
702  AliDebugStream(3) << "Creating axis " << axis.first << " with nBins " << axisLimits.at(0) << ", min: " << axisLimits.at(1) << ", max: " << axisLimits.at(2) << "\n";
703  binning.emplace_back(new TAxis(axisLimits.at(0), axisLimits.at(1), axisLimits.at(2)));
704  }
705 
706  // "s" ensures that Sumw2() is called
707  // The explicit const_cast is required
708  THnSparse * hist = fHistManager.CreateTHnSparse(name.c_str(), title.c_str(), binning.size(), const_cast<const TAxis **>(binning.data()), "s");
709  // Set the axis titles
710  int axisNumber = 0;
711  for (auto axis = sparseAxes.begin(); axis != sparseAxes.end(); axis++) {
712  AliDebugStream(5) << "ResponseMatrix: Add axis " << axis->first << " to sparse\n";
713  hist->GetAxis(axisNumber)->SetTitle(axis->first.c_str());
714  axisNumber++;
715  }
716 
717  // Define mapping from value name to value (including jet information)
718  for (unsigned int iJet = 1; iJet < 3; iJet++)
719  {
720  // pt
721  // ex: p_{T,1}
722  fResponseMatrixFillMap.insert(std::make_pair(std::string("p_{T,") + std::to_string(iJet) + "}", std::make_pair(iJet, &ResponseMatrixFillWrapper::fPt)));
723  // Area
724  // ex: A_{jet,1}
725  fResponseMatrixFillMap.insert(std::make_pair(std::string("A_{jet,") + std::to_string(iJet) + "}", std::make_pair(iJet, &ResponseMatrixFillWrapper::fArea)));
726  // EP angle
727  // ex: #theta_{jet,1}^{EP}
728  fResponseMatrixFillMap.insert(std::make_pair(std::string("#theta_{jet,") + std::to_string(iJet) + "}^{EP}", std::make_pair(iJet, &ResponseMatrixFillWrapper::fRelativeEPAngle)));
729  // Leading hadron
730  // ex: p_{T,particle,1}^{leading} (GeV/c)
731  fResponseMatrixFillMap.insert(std::make_pair(std::string("p_{T,particle,") + std::to_string(iJet) + "}^{leading} (GeV/c)", std::make_pair(iJet, &ResponseMatrixFillWrapper::fLeadingHadronPt)));
732  }
733  // Distance from one jet to another
734  fResponseMatrixFillMap.insert(std::make_pair("distance", std::make_pair(1, &ResponseMatrixFillWrapper::fDistance)) );
735  // Centrality
736  fResponseMatrixFillMap.insert(std::make_pair("centrality", std::make_pair(1, &ResponseMatrixFillWrapper::fCentrality)) );
737 
738  // Shared momentum fraction
739  name = "fHistFractionSharedPt";
740  title = "Fraction of p_{T} shared between matched jets";
741  // Need to include the bin from 1-1.01 to ensure that jets which shared all of their momentum
742  // due not end up in the overflow bin!
743  fHistManager.CreateTH1(name.c_str(), title.c_str(), 101, 0, 1.01);
744 }
745 
750 {
751  // Only fill the embedding qa plots if:
752  // - We are using the embedding helper
753  // - The class has been initialized
754  if (fEmbeddingQA.IsInitialized()) {
756  }
757 
758  // QA
759  if (fCreateQAHists) {
760  QAHists();
761  }
762 
763  // Jet matching
764  if (fPerformJetMatching) {
765  // Setup
766  AliDebugStream(1) << "Performing jet matching\n";
767  // We don't perform any additional jet matching initialization because we will restrict
768  // the jets accepted via the jet acceptance cuts (ie. EMCal fiducial cuts)
769  // Retrieve the releveant jet collections
770  AliJetContainer * jetsHybrid = GetJetContainer("hybridLevelJets");
771  AliJetContainer * jetsDetLevel = GetJetContainer("detLevelJets");
772  AliJetContainer * jetsPartLevel = GetJetContainer("partLevelJets");
773  // Validation
774  if (!jetsHybrid) {
775  AliErrorStream() << "Could not retrieve hybrid jet collection.\n";
776  return kFALSE;
777  }
778  if (!jetsDetLevel) {
779  AliErrorStream() << "Could not retrieve det level jet collection.\n";
780  return kFALSE;
781  }
782  if (!jetsPartLevel) {
783  AliErrorStream() << "Could not retrieve part level jet collection.\n";
784  return kFALSE;
785  }
786 
787  // Now, begin the actual matching.
788  // Hybrid <-> det first
789  AliDebugStream(2) << "Matching hybrid to detector level jets.\n";
790  // First, we reset the tagging
791  ResetMatching(*jetsHybrid);
792  ResetMatching(*jetsDetLevel);
793  // Next, we perform the matching
794  PerformGeometricalJetMatching(*jetsHybrid, *jetsDetLevel, fMaxJetMatchingDistance);
795  // Now, begin the next matching stage
796  // det <-> particle
797  AliDebugStream(2) << "Matching detector level to particle level jets.\n";
798  // First, we reset the tagging. We need to reset the det matching again to ensure
799  // that it doesn't accidentally keep some latent matches to the hybrid jets.
800  ResetMatching(*jetsDetLevel);
801  ResetMatching(*jetsPartLevel);
802  // Next, we perform the matching
803  PerformGeometricalJetMatching(*jetsDetLevel, *jetsPartLevel, fMaxJetMatchingDistance);
804 
805  // Fill QA hists.
806  FillJetMatchingQA(*jetsHybrid, *jetsDetLevel, "hybridToDet");
807  FillJetMatchingQA(*jetsDetLevel, *jetsPartLevel, "detToPart");
808  }
809 
810  // Response matrix
811  if (fCreateResponseMatrix) {
812  ResponseMatrix();
813  }
814 
815  return kTRUE;
816 }
817 
822 {
823  // Continue on to filling the histograms
824  FillQAHists();
825 }
826 
834 void AliAnalysisTaskEmcalJetHPerformance::FillCellQAHists(const std::string & prefix, AliVCaloCells * cells)
835 {
836  AliDebugStream(4) << "Storing cells with prefix \"" << prefix << "\". N cells: " << cells->GetNumberOfCells() << "\n";
837  short absId = -1;
838  double eCell = 0;
839  double tCell = 0;
840  double eFrac = 0;
841  int mcLabel = -1;
842 
843  std::string energyName = prefix + "/fHistCellEnergy";
844  std::string timeName = prefix + "/fHistCellTime";
845  std::string idName = prefix + "/fHistCellID";
846  bool embeddedCellWithLateCellTime = false;
847  bool fillingEmbeddedCells = (prefix.find("embedding") != std::string::npos);
848  for (unsigned int iCell = 0; iCell < cells->GetNumberOfCells(); iCell++) {
849  cells->GetCell(iCell, absId, eCell, tCell, mcLabel, eFrac);
850 
851  AliDebugStream(5) << "Cell " << iCell << ": absId: " << absId << ", E: " << eCell << ", t: " << tCell
852  << ", mcLabel: " << mcLabel << ", eFrac: " << eFrac << "\n";
853  fHistManager.FillTH1(energyName.c_str(), eCell);
854  fHistManager.FillTH1(timeName.c_str(), tCell);
855  fHistManager.FillTH1(idName.c_str(), absId);
856 
857  // We will record the event selection if the time is less than -400 ns
858  // This corresponds to a doubly corrected embedded event, which shouldn't be possible, and therefore
859  // indicates that something has gone awry
860  // NOTE: We must also require that the time is greater than -1 because apparently some uncalibrated cells
861  // will report their time as -1. We don't want to include those cells.
862  if (tCell < -400e-9 && tCell > -1 && fillingEmbeddedCells) {
863  embeddedCellWithLateCellTime = true;
864  }
865  }
866 
867  // If we have one embedded cell with late cell time, then we want to fill out the QA to
868  // help identify the event.
869  std::string embeddedEventUsed = prefix + "/fHistEmbeddedEventUsed";
870  std::string embeddedInteranlEventSelection = prefix + "/fHistInternalEventSelection";
871  if (embeddedCellWithLateCellTime)
872  {
873  auto embeddingInstance = AliAnalysisTaskEmcalEmbeddingHelper::GetInstance();
874  if (embeddingInstance) {
875  fHistManager.FillTH1(embeddedEventUsed.c_str(),
876  static_cast<double>(fPreviousEmbeddedEventSelected));
877 
878  // Determine the physics selection. This isn't quite a perfect way to store it, as it mingles the
879  // selections between different events. But it is simple, which will let us investigate quickly.
880  // Plus, it's a reasonable bet that the event selection when be the same when it goes wrong.
881  std::bitset<sizeof(UInt_t) * 8> testBits = fPreviousEventTrigger;
882  for (unsigned int i = 0; i < 32; i++) {
883  if (testBits.test(i)) {
884  fHistManager.FillTH2(embeddedInteranlEventSelection.c_str(),
885  static_cast<double>(fPreviousEmbeddedEventSelected), i);
886  }
887  }
888  }
889  }
890 
891  std::string nCellsName = prefix + "/fHistNCells";
892  fHistManager.FillTH1(nCellsName.c_str(), cells->GetNumberOfCells());
893 }
894 
899 {
900  // Fill standard cell QA
901  std::string prefix = "QA/";
902  prefix += fCaloCellsName.Data();
903  FillCellQAHists(prefix, fCaloCells);
904 
905  // Fill embedded cell QA it if's available.
906  auto embeddingInstance = AliAnalysisTaskEmcalEmbeddingHelper::GetInstance();
907  if (embeddingInstance) {
908  auto embeddedCells = dynamic_cast<AliVCaloCells*>(
909  embeddingInstance->GetExternalEvent()->FindListObject(fEmbeddedCellsName.c_str()));
910  if (embeddedCells) {
911  prefix = "QA/embedding/";
912  prefix += fEmbeddedCellsName;
913  FillCellQAHists(prefix, embeddedCells);
914  }
915  }
916 }
917 
922 {
923  if (fCreateCellQAHists) {
924  FillCellQAHists();
925  }
926 
927  // Tracks
928  AliTrackContainer * trackCont = 0;
929  TIter nextTrackColl(&fParticleCollArray);
930  while ((trackCont = static_cast<AliTrackContainer*>(nextTrackColl()))) {
931  for (auto track : trackCont->accepted())
932  {
933  fHistManager.FillTH3(TString::Format("QA/%s/fHistTrackPtEtaPhi", trackCont->GetName()), track->Pt(),
934  track->Eta(), track->Phi());
935  fHistManager.FillTH3(TString::Format("QA/%s/fHistTrackPtEtaPhiEfficiencyCorrected", trackCont->GetName()),
936  track->Pt(), track->Eta(), track->Phi(),
937  DetermineTrackingEfficiency(track->Pt(), track->Eta()));
938  }
939  }
940 
941  // Clusters
942  AliClusterContainer* clusCont = 0;
943  TIter nextClusColl(&fClusterCollArray);
944  AliVCluster * cluster = 0;
945  while ((clusCont = static_cast<AliClusterContainer*>(nextClusColl()))) {
946  for (auto clusIter : clusCont->accepted_momentum())
947  {
948  AliTLorentzVector c = clusIter.first;
949  cluster = clusIter.second;
950  // Intentionally plotting against raw energy
951  fHistManager.FillTH2(TString::Format("QA/%s/fHistClusterEnergyVsTime", clusCont->GetName()), cluster->E(), cluster->GetTOF());
952  // Hadronically corrected energy (which is usually what we're using)
953  fHistManager.FillTH1(TString::Format("QA/%s/fHistClusterHadCorrEnergy", clusCont->GetName()),
954  cluster->GetHadCorrEnergy());
955  fHistManager.FillTH2(TString::Format("QA/%s/fHistClusterEtaPhi", clusCont->GetName()), c.Eta(),
956  c.Phi_0_2pi());
957  }
958  }
959 
960  // Jets
961  AliJetContainer * jetCont = 0;
962  double jetPt = 0;
963  TIter nextJetColl(&fJetCollArray);
964  while ((jetCont = static_cast<AliJetContainer*>(nextJetColl()))) {
965  double rhoVal = jetCont->GetRhoVal();
966  for (auto jet : jetCont->accepted())
967  {
968  jetPt = AliAnalysisTaskEmcalJetHUtils::GetJetPt(jet, rhoVal);
969  fHistManager.FillTH1(TString::Format("QA/%s/fHistJetPt", jetCont->GetName()), jetPt);
970  }
971  }
972 
973  // Event plane resolution
974  // Only attempt to fill the histogram if the corrected Qn vectors are available.
975  if (fFlowQnVectorManager) {
976  // Retrieve the Qn vectors
977  const AliQnCorrectionsQnVector * vzeroQnVector = fFlowQnVectorManager->GetDetectorQnVector("VZEROQoverM");
978  const AliQnCorrectionsQnVector * tpcPosQnVector = fFlowQnVectorManager->GetDetectorQnVector("TPCPosEtaQoverM");
979  const AliQnCorrectionsQnVector * tpcNegQnVector = fFlowQnVectorManager->GetDetectorQnVector("TPCNegEtaQoverM");
980  if (vzeroQnVector == nullptr || tpcPosQnVector == nullptr || tpcNegQnVector == nullptr) {
981  AliWarningStream() << "Q vector unavailable. Skipping.\n";
982  }
983  else {
984  // We need the centrality bin, but the bin provided by the base class isn't that fine grained
985  // (and may cause problems elsewhere if it is). So we just calculate the value here.
986  std::string name = "QA/eventPlaneRes/%s/R%i";
987  // The resolution is calculated by changing the leading term, but it is all with respect
988  // to the same harmonic. We measure with respect to the second order harmonic, so we retrieve
989  // that event plane.
990  double vzero = vzeroQnVector->EventPlane(2);
991  double tpcPos = tpcPosQnVector->EventPlane(2);
992  double tpcNeg = tpcNegQnVector->EventPlane(2);
993  for (unsigned int R = 2; R < 9; R++) {
994  fHistManager.FillProfile(TString::Format(name.c_str(), "VZERO", R), fCent,
995  std::cos(R * (tpcPos - tpcNeg)));
996  fHistManager.FillProfile(TString::Format(name.c_str(), "TPC_Pos", R), fCent,
997  std::cos(R * (vzero - tpcNeg)));
998  fHistManager.FillProfile(TString::Format(name.c_str(), "TPC_Neg", R), fCent,
999  std::cos(R * (vzero - tpcPos)));
1000  }
1001  }
1002  }
1003 
1004  // Update the previous event trigger to the current event trigger so that it is available next event.
1005  // This is stored for keeping track of when embedded events are double corrected.
1006  // This must be updated after filling the relevant hists above!
1008  ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
1009  auto embeddingInstance = AliAnalysisTaskEmcalEmbeddingHelper::GetInstance();
1010  if (embeddingInstance) {
1011  fPreviousEmbeddedEventSelected = embeddingInstance->EmbeddedEventUsed();
1012  }
1013 }
1014 
1021 {
1022  for(auto j : c.all()){
1023  j->ResetMatching();
1024  }
1025 }
1026 
1050  AliJetContainer& contTag, double maxDist) const
1051 {
1052  // Setup
1053  const Int_t kNacceptedBase = contBase.GetNAcceptedJets(), kNacceptedTag = contTag.GetNAcceptedJets();
1054  if (!(kNacceptedBase && kNacceptedTag)) {
1055  return false;
1056  }
1057 
1058  // Build up vectors of jet pointers to use when assigning the closest jets.
1059  // The storages are needed later for applying the tagging, in order to avoid multiple occurrence of jet selection
1060  std::vector<AliEmcalJet*> jetsBase(kNacceptedBase), jetsTag(kNacceptedTag);
1061 
1062  int countBase(0), countTag(0);
1063  for (auto jb : contBase.accepted()) {
1064  jetsBase[countBase] = jb;
1065  countBase++;
1066  }
1067  for (auto jt : contTag.accepted()) {
1068  jetsTag[countTag] = jt;
1069  countTag++;
1070  }
1071 
1072  TArrayI faMatchIndexTag(kNacceptedBase), faMatchIndexBase(kNacceptedTag);
1073  faMatchIndexBase.Reset(-1);
1074  faMatchIndexTag.Reset(-1);
1075 
1076  // find the closest distance to the base jet
1077  countBase = 0;
1078  for (auto jet1 : contBase.accepted()) {
1079  double distance = maxDist;
1080 
1081  // Loop over all accepted jets and brute force search for the closest jet.
1082  // NOTE: current_index() returns the jet index in the underlying array, not
1083  // the index within the accepted jets that are returned.
1084  int contTagAcceptedIndex = 0;
1085  for (auto jet2 : contTag.accepted()) {
1086  double dR = jet1->DeltaR(jet2);
1087  if (dR < distance && dR < maxDist) {
1088  faMatchIndexTag[countBase] = contTagAcceptedIndex;
1089  distance = dR;
1090  }
1091  contTagAcceptedIndex++;
1092  }
1093 
1094  // Let us know whether a match was found successfully.
1095  if (faMatchIndexTag[countBase] >= 0 && distance < maxDist) {
1096  AliDebugStream(1) << "Found closest tag jet for " << countBase << " with match index "
1097  << faMatchIndexTag[countBase] << " and distance " << distance << "\n";
1098  } else {
1099  AliDebugStream(1) << "Not found closest tag jet for " << countBase << ".\n";
1100  }
1101 
1102  countBase++;
1103  }
1104 
1105  // other way around
1106  countTag = 0;
1107  for (auto jet1 : contTag.accepted()) {
1108  double distance = maxDist;
1109 
1110  // Loop over all accepted jets and brute force search for the closest jet.
1111  // NOTE: current_index() returns the jet index in the underlying array, not
1112  // the index within the accepted jets that are returned.
1113  int contBaseAcceptedIndex = 0;
1114  for (auto jet2 : contBase.accepted()) {
1115  double dR = jet1->DeltaR(jet2);
1116  if (dR < distance && dR < maxDist) {
1117  faMatchIndexBase[countTag] = contBaseAcceptedIndex;
1118  distance = dR;
1119  }
1120  contBaseAcceptedIndex++;
1121  }
1122 
1123  // Let us know whether a match was found successfully.
1124  if (faMatchIndexBase[countTag] >= 0 && distance < maxDist) {
1125  AliDebugStream(1) << "Found closest base jet for " << countTag << " with match index "
1126  << faMatchIndexBase[countTag] << " and distance " << distance << "\n";
1127  } else {
1128  AliDebugStream(1) << "Not found closest base jet for " << countTag << ".\n";
1129  }
1130 
1131  countTag++;
1132  }
1133 
1134  // check for "true" correlations
1135  // these are pairs where the base jet is the closest to the tag jet and vice versa
1136  // As the lists are linear a loop over the outer base jet is sufficient.
1137  AliDebugStream(1) << "Starting true jet loop: nbase(" << kNacceptedBase << "), ntag(" << kNacceptedTag << ")\n";
1138  for (int ibase = 0; ibase < kNacceptedBase; ibase++) {
1139  AliDebugStream(2) << "base jet " << ibase << ": match index in tag jet container " << faMatchIndexTag[ibase]
1140  << "\n";
1141  if (faMatchIndexTag[ibase] > -1) {
1142  AliDebugStream(2) << "tag jet " << faMatchIndexTag[ibase] << ": matched base jet "
1143  << faMatchIndexBase[faMatchIndexTag[ibase]] << "\n";
1144  }
1145  // We have a true correlation where each jet points to the other.
1146  if (faMatchIndexTag[ibase] > -1 && faMatchIndexBase[faMatchIndexTag[ibase]] == ibase) {
1147  AliDebugStream(2) << "found a true match \n";
1148  AliEmcalJet *jetBase = jetsBase[ibase], *jetTag = jetsTag[faMatchIndexTag[ibase]];
1149  // We have a valid pair of matched jets, so set the closest jet properties.
1150  if (jetBase && jetTag) {
1151  Double_t dR = jetBase->DeltaR(jetTag);
1152  jetBase->SetClosestJet(jetTag, dR);
1153  jetTag->SetClosestJet(jetBase, dR);
1154  }
1155  }
1156  }
1157  return true;
1158 }
1159 
1167  const std::string& prefix)
1168 {
1169  // Fill histograms.
1170  std::string name = "";
1171  std::string centBin = std::to_string(fCentBin);
1172  AliDebugStream(2) << "Filling matching hists with prefix: " << prefix << "\n";
1173  for (auto jet1 : contBase.accepted()) {
1174  // Setup
1175  Double_t ptJet1 = jet1->Pt() - contBase.GetRhoVal() * jet1->Area();
1176  // Record jet 1 only properties
1177  AliDebugStream(4) << "jet1: " << jet1->toString() << "\n";
1178  name = "jetMatching/" + prefix + "/fh2PtJet1VsLeadPtAllSel_" + centBin;
1179  fHistManager.FillTH2(name.c_str(), ptJet1, jet1->MaxTrackPt());
1180 
1181  // Retrieve jet 2
1182  AliEmcalJet * jet2 = jet1->ClosestJet();
1183  if (!jet2) { continue; }
1184  AliDebugStream(4) << "jet2: " << jet2->toString() << "\n";
1185  Double_t ptJet2 = jet2->Pt() - contTag.GetRhoVal() * jet2->Area();
1186  // This will retrieve the fraction of jet2's momentum in jet1.
1187  Double_t fraction = contBase.GetFractionSharedPt(jet1);
1188 
1189  name = "jetMatching/" + prefix + "/fh2PtJet2VsFraction_" + centBin;
1190  fHistManager.FillTH2(name.c_str(), ptJet2, fraction);
1191  AliDebugStream(5) << "Fraction: " << fraction << ", minimum: " << fMinFractionShared << "\n";
1192  if (fraction < fMinFractionShared) {
1193  continue;
1194  }
1195  name = "jetMatching/" + prefix + "/fh2PtJet1VsLeadPtTagged_" + centBin;
1196  fHistManager.FillTH2(name.c_str(), ptJet1,jet1->MaxTrackPt());
1197  name = "jetMatching/" + prefix + "/fh2PtJet1VsPtJet2_" + centBin;
1198  fHistManager.FillTH2(name.c_str(), ptJet1, ptJet2);
1199  if (ptJet2 > 0.) {
1200  name = "jetMatching/" + prefix + "/fh2PtJet2VsRelPt_" + centBin;
1201  fHistManager.FillTH2(name.c_str(), ptJet2, (ptJet1 - ptJet2) / ptJet2);
1202  }
1203 
1204  // Recall that the arguments are backwards of the expectation.
1205  Double_t dPhi = DeltaPhi(jet2->Phi(), jet1->Phi());
1206  if (dPhi > TMath::Pi()) {
1207  dPhi -= TMath::TwoPi();
1208  }
1209  if (dPhi < (-1. * TMath::Pi())) {
1210  dPhi += TMath::TwoPi();
1211  }
1212 
1213  name = "jetMatching/" + prefix + "/fh3PtJet1VsDeltaEtaDeltaPhi_" + centBin;
1214  fHistManager.FillTH3(name.c_str(), ptJet1, jet1->Eta() - jet2->Eta(), dPhi);
1215  name = "jetMatching/" + prefix + "/fh2PtJet1VsDeltaR_" + centBin;
1216  fHistManager.FillTH2(name.c_str(), ptJet1, jet1->DeltaR(jet2));
1217  }
1218 
1219  // Number of accepted jets
1220  name = "jetMatching/" + prefix + "/fNAccJets";
1221  fHistManager.FillTH1(name.c_str(), contBase.GetNAcceptedJets());
1222 }
1223 
1228 {
1229  AliJetContainer * jetsHybrid = GetJetContainer("hybridLevelJets");
1230  AliJetContainer * jetsDetLevel = GetJetContainer("detLevelJets");
1231  AliJetContainer * jetsPartLevel = GetJetContainer("partLevelJets");
1232  // NOTE: Defaults to 0 if rho was not specified.
1233  double rhoHybrid = jetsHybrid->GetRhoVal();
1234  if (!jetsHybrid) {
1235  AliErrorStream() << "Could not retrieve hybrid jet collection.\n";
1236  return;
1237  }
1238  if (!jetsDetLevel) {
1239  AliErrorStream() << "Could not retrieve det level jet collection.\n";
1240  return;
1241  }
1242  if (fResponseFromThreeJetCollections && !jetsPartLevel) {
1243  AliErrorStream() << "Could not retrieve part level jet collection.\n";
1244  return;
1245  }
1246 
1247  // Handle matching of jets.
1248  for (auto jet1 : jetsHybrid->accepted())
1249  {
1250  AliDebugStream(4) << "jet1: " << jet1->toString() << "\n";
1251  AliDebugStream(4) << "jet1 address: " << jet1 << "\n";
1252 
1253  // Get jet the det level jet from the hybrid jet
1254  AliEmcalJet * jet2 = jet1->ClosestJet();
1255  if(!jet2) continue;
1256 
1257  AliDebugStream(4) << "jet2: " << jet2->toString() << "\n";
1258  AliDebugStream(4) << "jet2 address: " << jet2 << "\n";
1259 
1260  // Check shared fraction
1261  double sharedFraction = jetsHybrid->GetFractionSharedPt(jet1);
1262  fHistManager.FillTH1("fHistFractionSharedPt", sharedFraction);
1263  if (sharedFraction < fMinFractionShared) {
1264  AliDebugStream(4) << "Rejecting jet due to momentum fraction of " << sharedFraction << ", smaller than the minimum.\n";
1265  continue;
1266  }
1267  else {
1268  AliDebugStream(4) << "Jet passed momentum fraction cut with value of " << sharedFraction << "\n";
1269  }
1270 
1271  // NOTE: We apply no explicit event selection to jet 2 because the matching in the tagger
1272  // only matches jets which are accepted.
1273 
1274  // Determine the jet to use to fill the response
1275  // The jet that is passed may be the embedded detector level (for two stage matching), or it may
1276  // be the embedded particle level (for direct matching).
1277  AliEmcalJet * jetToPass = 0;
1279  // Retrieve the MC level jet
1280  AliEmcalJet * jet3 = jet2->ClosestJet();
1281  UInt_t rejectionReason = 0;
1282  if (!jetsPartLevel->AcceptJet(jet3, rejectionReason)) {
1283  // NOTE: This shouldn't ever happen because the tagger applies acceptance
1284  // cuts when matching. However, we keep the check here for good measure.
1285  // NOTE: Could store rejection reasons if needed with below:
1286  //fHistRejectionReason2->Fill(jets2->GetRejectionReasonBitPosition(rejectionReason), jet2->Pt());
1287  continue;
1288  }
1289 
1290  AliDebugStream(4) << "jet3: " << jet3->toString() << "\n";
1291  AliDebugStream(4) << "jet3 address: " << jet3 << "\n";
1292 
1293  // Use for the response
1294  AliDebugStream(4) << "Using part level jet for response (ie. two stage matching)\n";
1295  jetToPass = jet3;
1296  }
1297  else {
1298  // Use for the response
1299  AliDebugStream(4) << "Using one stage matching for response\n";
1300  jetToPass = jet2;
1301  }
1302 
1303  // Fill response
1304  FillResponseMatrix(jet1, jetToPass, rhoHybrid);
1305  }
1306 
1307 }
1308 
1313 {
1314  if (!jet1 || !jet2) {
1315  AliErrorStream() << "Null jet passed to fill response matrix";
1316  }
1317 
1318  AliDebugStream(3) << "About to create ResponseMatrixFillWrappers\n";
1319  AliDebugStream(4) << "jet1: " << jet1->toString() << "\n";
1320  AliDebugStream(4) << "jet2: " << jet2->toString() << "\n";
1321  // Create map from jetNumber to jet and initialize the objects
1322  std::map<unsigned int, ResponseMatrixFillWrapper> jetNumberToJet = {
1323  std::make_pair(1, CreateResponseMatrixFillWrapper(jet1, jet1Rho)),
1324  // We would never want to rho subtract the second jet, regardless of whether it's external
1325  // detector level or particle level.
1326  std::make_pair(2, CreateResponseMatrixFillWrapper(jet2, 0))
1327  };
1328 
1329  // Fill histograms
1330  std::string histName = "response/fHistResponseMatrix";
1331  std::vector<double> values;
1332  THnSparse * response = static_cast<THnSparse*>(fHistManager.FindObject(histName.c_str()));
1333  AliDebugStream(3) << "About to fill response matrix values\n";
1334  AliDebugStream(4) << "jet1: " << jet1->toString() << "\n";
1335  AliDebugStream(4) << "jet2: " << jet2->toString() << "\n";
1336  for (unsigned int i = 0; i < response->GetNdimensions(); i++) {
1337  std::string title = response->GetAxis(i)->GetTitle();
1338 
1339  // Retrieve pair of jet and pointer to extract the fill value
1340  auto jetPair = fResponseMatrixFillMap.find(title);
1341  if (jetPair != fResponseMatrixFillMap.end()) {
1342  auto wrapper = jetNumberToJet.at(jetPair->second.first);
1343  auto member = jetPair->second.second;
1344  AliDebugStream(4) << "Filling value " << wrapper.*member << " into axis " << title << "\n";
1345  values.emplace_back(wrapper.*member);
1346  }
1347  else {
1348  AliWarningStream() << "Unable to fill dimension " << title << "!\n";
1349  }
1350  }
1351 
1352  fHistManager.FillTHnSparse(histName.c_str(), values.data());
1353 }
1354 
1359 {
1360  ResponseMatrixFillWrapper wrapper;
1361  if (!jet) {
1362  AliErrorStream() << "Must pass valid jet to create object.\n";
1363  return wrapper;
1364  }
1365  wrapper.fPt = AliAnalysisTaskEmcalJetHUtils::GetJetPt(jet, rho);
1366  wrapper.fArea = jet->Area();
1367  wrapper.fPhi = jet->Phi();
1368  wrapper.fDistance = jet->ClosestJetDistance();
1369  wrapper.fCentrality = fCent;
1372 
1373  return wrapper;
1374 }
1375 
1382 {
1383  // Get the pointer to the existing analysis manager via the static access method.
1384  //==============================================================================
1385  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1386  if (!mgr)
1387  {
1388  AliErrorClass("No analysis manager to connect to.");
1389  return nullptr;
1390  }
1391 
1392  // Setup task name
1393  std::string taskName = "AliAnalysisTaskEmcalJetHPerformance";
1394  std::string suffixName(suffix);
1395  if (suffixName != "") {
1396  taskName += "_";
1397  taskName += suffixName;
1398  }
1399 
1400  // Create task and configure as desired.
1402  // Set a few general default.
1403  task->SetNCentBins(5);
1404  // Configuration is via YAML.
1405  mgr->AddTask(task);
1406 
1407  // Create containers for input/output
1408  mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer() );
1409  AliAnalysisDataContainer* outputContainer =
1410  mgr->CreateContainer(task->GetName(), TList::Class(), AliAnalysisManager::kOutputContainer,
1411  Form("%s", AliAnalysisManager::GetCommonFileName()));
1412  mgr->ConnectOutput(task, 1, outputContainer);
1413 
1414  return task;
1415 }
1416 
1423 {
1424  std::stringstream tempSS;
1425  tempSS << std::boolalpha;
1426  tempSS << "Particle collections:\n";
1427  TIter nextParticleCont(&fParticleCollArray);
1428  AliParticleContainer * particleCont;
1429  while ((particleCont = static_cast<AliParticleContainer *>(nextParticleCont()))) {
1430  tempSS << "\t" << particleCont->GetName() << ": " << particleCont->GetArrayName() << "\n";
1431  }
1432  tempSS << "Cluster collections:\n";
1433  TIter nextClusterCont(&fClusterCollArray);
1434  AliClusterContainer * clusterCont;
1435  while ((clusterCont = static_cast<AliClusterContainer *>(nextClusterCont()))) {
1436  tempSS << "\t" << clusterCont->GetName() << ": " << clusterCont->GetArrayName() << "\n";
1437  }
1438  tempSS << "Jet collections:\n";
1439  TIter nextJetCont(&fJetCollArray);
1440  AliJetContainer * jetCont;
1441  while ((jetCont = static_cast<AliJetContainer *>(nextJetCont()))) {
1442  tempSS << "\t" << jetCont->GetName() << ": " << jetCont->GetArrayName() << "\n";
1443  }
1444  // AliEventCuts
1445  tempSS << "AliEventCuts\n";
1446  tempSS << "\tEnabled: " << !fUseBuiltinEventSelection << "\n";
1447  // Efficiency
1448  tempSS << "Efficiency\n";
1449  tempSS << "\tSingle track efficiency identifier: " << fEfficiencyPeriodIdentifier << "\n";
1450  // QA
1451  tempSS << "QA Hists:\n";
1452  tempSS << "\tEnabled: " << fCreateQAHists << "\n";
1453  tempSS << "\tCreate cell QA hists: " << fCreateCellQAHists << "\n";
1454  // Use whether the pointer as null as a proxy. It's not ideal because it's not fully initialized
1455  // until after UserCreateOutputObjects(). But it's good enough for these purposes.
1456  tempSS << "\tCalculate event plane resolution (proxy of whether it's enabled - it may not be accurate): " << (fFlowQnVectorManager != nullptr) << "\n";
1457  // Jet matching
1458  tempSS << "Jet matching:\n";
1459  tempSS << "\tEnabled: " << fPerformJetMatching << "\n";
1460  tempSS << "\tMax matching distance: " << fMaxJetMatchingDistance << "\n";
1461  // Response matrix
1462  tempSS << "Response matrix:\n";
1463  tempSS << "\tEnabled: " << fCreateResponseMatrix << "\n";
1464  tempSS << "\tConstruct response from 3 jet collections: " << fResponseFromThreeJetCollections << "\n";
1465  tempSS << "\tMin fraction shared pt: " << fMinFractionShared << "\n";
1466  tempSS << "\tJet leading hadron bias type: " << fLeadingHadronBiasType << "\n";
1467  tempSS << "\tResponse matrix fill map: \n";
1468  for (auto el : fResponseMatrixFillMap) {
1469  tempSS << "\t\tProperty " << el.first << " applied to jet " << el.second.first << "\n";
1470  }
1471 
1472  return tempSS.str();
1473 }
1474 
1481 std::ostream & AliAnalysisTaskEmcalJetHPerformance::Print(std::ostream & in) const {
1482  in << toString();
1483  return in;
1484 }
1485 
1493 {
1494  Printf("%s", toString().c_str());
1495 }
1496 
1497 } /* namespace EMCALJetTasks */
1498 } /* namespace PWGJE */
1499 
1507 std::ostream & operator<<(std::ostream & in, const PWGJE::EMCALJetTasks::AliAnalysisTaskEmcalJetHPerformance & myTask)
1508 {
1509  std::ostream & result = myTask.Print(in);
1510  return result;
1511 }
1512 
1519 {
1520  using std::swap;
1521 
1522  // Same ordering as in the constructors (for consistency)
1523  swap(first.fYAMLConfig, second.fYAMLConfig);
1525  //swap(first.fHistManager, second.fHistManager); // Skip here, because the THistManager copy constructor is private.
1526  //swap(first.fEmbeddingQA, second.fEmbeddingQA); // Skip here, because the THistManager copy constructor is private.
1527  swap(first.fCreateQAHists, second.fCreateQAHists);
1540 }
static void ConfigureTrackContainersFromYAMLConfig(std::vector< std::string > baseNameWithContainer, AliTrackContainer *trackCont, PWG::Tools::AliYAMLConfiguration &yamlConfig, std::string taskName)
void FillResponseMatrix(AliEmcalJet *jet1, AliEmcalJet *jet2, const double jet1Rho)
Double_t Area() const
Definition: AliEmcalJet.h:130
TObjArray fClusterCollArray
cluster collection array
Double_t GetRhoVal() const
double Double_t
Definition: External.C:58
static double DetermineTrackingEfficiency(const double trackPt, const double trackEta, const int centralityBin, const EEfficiencyPeriodIdentifier_t efficiencyPeriodIdentifier, const std::string &taskName)
bool fResponseFromThreeJetCollections
If true, the det level jets in collection 2 are only an intermediate step. They are used to get part ...
static Double_t DeltaPhi(Double_t phia, Double_t phib, Double_t rMin=-TMath::Pi()/2, Double_t rMax=3 *TMath::Pi()/2)
const char * title
Definition: MakeQAPdf.C:27
AliEmcalJet * ClosestJet() const
Definition: AliEmcalJet.h:341
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)
void SetName(const char *n)
Set the name of the class of the objets inside the underlying array.
Double_t ClosestJetDistance() const
Definition: AliEmcalJet.h:342
ResponseMatrixFillWrapper CreateResponseMatrixFillWrapper(AliEmcalJet *jet, const double rho) const
void AdoptParticleContainer(AliParticleContainer *cont)
Double_t Eta() const
Definition: AliEmcalJet.h:121
const AliClusterIterableMomentumContainer accepted_momentum() const
void SetLeadingHadronType(Int_t t)
const AliTrackIterableContainer accepted() const
Double_t Phi() const
Definition: AliEmcalJet.h:117
Container with name, TClonesArray and cuts for particles.
void FillTH2(const char *hname, double x, double y, double weight=1., Option_t *opt="")
Fill a 2D histogram within the container.
Declaration of class AliTLorentzVector.
double fCentrality
! Centrality of the given jet (an event level property, but useful to here available here) ...
Double_t fEPV0
!event plane V0
void FillTH3(const char *hname, double x, double y, double z, double weight=1., Option_t *opt="")
Fill a 3D histogram within the container.
Declaration of class AliAnalysisTaskEmcalEmbeddingHelper.
AliAnalysisTaskEmcalJetHUtils::EEfficiencyPeriodIdentifier_t fEfficiencyPeriodIdentifier
Identifies the period for determining the efficiency correction to apply.
virtual Bool_t AcceptJet(Int_t i, UInt_t &rejectionReason) const
TCanvas * c
Definition: TestFitELoss.C:172
AliJetContainer * AddJetContainer(const char *n, TString defaultCutType, Float_t jetRadius=0.4)
AliAnalysisTaskEmcalJetHPerformance & operator=(AliAnalysisTaskEmcalJetHPerformance other)
Int_t fCentBin
!event centrality bin
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 SetPercAreaCut(Float_t p)
void SetObject(TObject *const o, const char *group="/")
Set a new group into the container into the parent group.
PWG::Tools::AliYAMLConfiguration fYAMLConfig
YAML configuration file.
bool fCreateResponseMatrix
If true, create a response matrix with the available jet collections.
Bool_t fUseBuiltinEventSelection
Use builtin event selection of the AliAnalysisTaskEmcal instead of AliEventCuts.
bool DoesConfigurationExist(const int i) const
AliClusterContainer * AddClusterContainer(const char *n)
Create new cluster container and attach it to the task.
Container for particles within the EMCAL framework.
const TString & GetClassName() const
TObjArray fParticleCollArray
particle/track collection array
bool AddQAPlotsToList(TList *list)
std::string fEmbeddedCellsName
Set the embedded cells collection name.
void SetRhoName(const char *n)
AliAnalysisTaskEmcalJetHUtils::ELeadingHadronBiasType_t fLeadingHadronBiasType
Leading hadron in jet bias type (either charged, neutral, or both)
static void ConfigureEventCuts(AliEventCuts &eventCuts, PWG::Tools::AliYAMLConfiguration &yamlConfig, const UInt_t offlineTriggerMask, const std::string &baseName, const std::string &taskName)
void SetMinTrackPt(Float_t b)
void FillTHnSparse(const char *name, const double *x, double weight=1., Option_t *opt="")
bool IsInitialized() const
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.
static UInt_t DeterminePhysicsSelectionFromYAML(const std::vector< std::string > &selections)
TObject * FindObject(const char *name) const
Find an object inside the container.
int Int_t
Definition: External.C:63
void CreateTProfile(const char *name, const char *title, int nbinsX, double xmin, double xmax, Option_t *opt="")
Create a new TProfile within 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
Double_t GetMinPt() const
UInt_t fPreviousEventTrigger
Physics selection (offline trigger) of the previous event for determine why a small number of embedde...
Double_t Phi_0_2pi() const
bool PerformGeometricalJetMatching(AliJetContainer &contBase, AliJetContainer &contTag, double maxDist) const
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.
static void ConfigureClusterContainersFromYAMLConfig(std::vector< std::string > baseNameWithContainer, AliClusterContainer *clusterCont, PWG::Tools::AliYAMLConfiguration &yamlConfig, std::string taskName)
Int_t fNcentBins
how many centrality bins
bool fCreateCellQAHists
If true, create the Cell QA histograms. It doesn&#39;t gracefully turn off when not configured like the c...
void FillProfile(const char *name, double x, double y, double weight=1.)
static double GetJetPt(const AliEmcalJet *jet, const double rho)
Double_t fCent
!event centrality
TString fCaloCellsName
name of calo cell collection
static AliParticleContainer * CreateParticleOrTrackContainer(const std::string &collectionName)
static const std::map< std::string, EEfficiencyPeriodIdentifier_t > fgkEfficiencyPeriodIdentifier
! Map from name to efficiency period identifier for use with the YAML config
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 void ConfigureEMCalContainersFromYAMLConfig(std::vector< std::string > baseName, std::string containerName, AliEmcalContainer *cont, PWG::Tools::AliYAMLConfiguration &yamlConfig, std::string taskName)
static double RelativeEPAngle(double jetAngle, double epAngle)
TObjArray fJetCollArray
jet collection array
Double_t DeltaR(const AliVParticle *part) const
virtual void SetNCentBins(Int_t n)
AliVCaloCells * fCaloCells
!cells
static double GetLeadingHadronPt(AliEmcalJet *jet, ELeadingHadronBiasType_t leadingHadronType)
Double_t Pt() const
Definition: AliEmcalJet.h:109
void FillJetMatchingQA(AliJetContainer &contBase, AliJetContainer &contTag, const std::string &prefix)
const TString & GetArrayName() const
bool fPreviousEmbeddedEventSelected
True if the previous embedded event was selected. Used to determine why a small number of embedded ev...
AliQnCorrectionsManager * fFlowQnVectorManager
! Qn corrections framework manager.
const char * GetName() const
static UInt_t DetermineJetAcceptanceFromYAML(const std::vector< std::string > &selections)
AliEventCuts fAliEventCuts
Event cuts (run2 defaults)
AliEmcalList * fOutput
!output list
void SetClosestJet(AliEmcalJet *j, Double_t d)
Definition: AliEmcalJet.h:337
void SetMakeGeneralHistograms(Bool_t g)
Enable general histograms.
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:51
static std::string DetermineUseDefaultName(InputObject_t objType)
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.
Container structure for EMCAL clusters.
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)
TH3 * CreateTH3(const char *name, const char *title, int nbinsx, double xmin, double xmax, int nbinsy, double ymin, double ymax, int nbinsz, double zmin, double zmax, Option_t *opt="")
Create a new TH2 within the container.
const AliJetIterableContainer all() const
static const AliAnalysisTaskEmcalEmbeddingHelper * GetInstance()