AliPhysics  64bcec2 (64bcec2)
AliAnalysisTaskEmcalEmbeddingHelper.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2016, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 #include <string>
17 #include <sstream>
18 #include <vector>
19 #include <algorithm>
20 #include <memory>
21 #include <fstream>
22 #include <iostream>
23 #include <bitset>
24 
25 #include <TFile.h>
26 #include <TMath.h>
27 #include <TRandom.h>
28 #include <TChain.h>
29 #include <TGrid.h>
30 #include <TGridResult.h>
31 #include <TSystem.h>
32 #include <TUUID.h>
33 #include <TKey.h>
34 #include <TProfile.h>
35 #include <TH1F.h>
36 #include <TRandom3.h>
37 #include <TList.h>
38 
39 #include <AliLog.h>
40 #include <AliAnalysisManager.h>
41 #include <AliVEvent.h>
42 #include <AliAODEvent.h>
43 #include <AliESDEvent.h>
44 #include <AliMCEvent.h>
45 #include <AliInputEventHandler.h>
46 #include <AliVHeader.h>
47 #include <AliAODMCHeader.h>
48 #include <AliGenPythiaEventHeader.h>
49 
50 #include "AliYAMLConfiguration.h"
51 #include "AliEmcalList.h"
52 #include "AliEmcalContainerUtils.h"
53 
55 
59 
64 {
65  if (!gGrid) {
66  AliInfoGeneralStream("AliAnalysisTaskEmcalEmbeddingHelper") << "Trying to connect to AliEn ...\n";
67  TGrid::Connect("alien://");
68  }
69  if (!gGrid) {
70  AliFatalGeneral("AliAnalysisTaskEmcalEmbeddingHelper", "Cannot access AliEn!");
71  }
72 }
73 
80 bool IsFileAccessible(std::string filename)
81 {
82  // Connect to AliEn if necessary
83  // Usually, `gGrid` will exist, so we won't need to waste time on find()
84  if (!gGrid && filename.find("alien://") != std::string::npos) {
86  }
87 
88  // AccessPathName() cannot handle the "#", so we need to strip it to check that the file exists.
89  if (filename.find(".zip#") != std::string::npos) {
90  std::size_t pos = filename.find_last_of("#");
91  filename.erase(pos);
92  }
93 
94  // AccessPathName() has an odd return value - false means that the file exists.
95  // NOTE: This is extremely inefficienct for TAlienSystem. It calls
96  // -> gapi_access() -> gapi_stat(). gapi_stat() calls "ls -l" on the basename directory,
97  // which can cause a load on AliEn.
98  bool res = gSystem->AccessPathName(filename.c_str());
99  // Normalize the result to true if file exists (needed because of the odd return value)
100  res = (res == false);
101  if (res == false) {
102  AliDebugGeneralStream("AliAnalysisTaskEmcalEmbeddingHelper", 4) << "File \"" << filename << "\" doees not exist!\n";
103  }
104 
105  return res;
106 }
107 
109 
115  fTriggerMask(0),
116  fMCRejectOutliers(false),
117  fPtHardJetPtRejectionFactor(4),
118  fZVertexCut(10),
119  fMaxVertexDist(999),
120  fInitializedConfiguration(false),
121  fInitializedNewFile(false),
122  fInitializedEmbedding(false),
123  fWrappedAroundTree(false),
124  fTreeName(""),
125  fNPtHardBins(1),
126  fPtHardBin(-1),
127  fRandomEventNumberAccess(kFALSE),
128  fRandomFileAccess(kTRUE),
129  fCreateHisto(true),
130  fYAMLConfig(),
131  fUseInternalEventSelection(false),
132  fUseManualInternalEventCuts(false),
133  fInternalEventCuts(),
134  fEmbeddedEventUsed(true),
135  fCentMin(-999),
136  fCentMax(-999),
137  fAutoConfigurePtHardBins(false),
138  fAutoConfigureBasePath(""),
139  fAutoConfigureTrainTypePath(""),
140  fAutoConfigureIdentifier(""),
141  fFilePattern(""),
142  fInputFilename(""),
143  fFileListFilename(""),
144  fFilenameIndex(-1),
145  fFilenames(),
146  fConfigurationPath(""),
147  fEmbeddedRunlist(),
148  fPythiaCrossSectionFilenames(),
149  fExternalFile(nullptr),
150  fChain(nullptr),
151  fCurrentEntry(0),
152  fLowerEntry(0),
153  fUpperEntry(0),
154  fOffset(0),
155  fMaxNumberOfFiles(0),
156  fFileNumber(0),
157  fHistManager(),
158  fOutput(nullptr),
159  fExternalEvent(nullptr),
160  fExternalHeader(nullptr),
161  fPythiaHeader(nullptr),
162  fPythiaTrials(0),
163  fPythiaTrialsFromFile(0),
164  fPythiaCrossSection(0.),
165  fPythiaCrossSectionFromFile(0.),
166  fPythiaPtHard(0.)
167 {
168  if (fgInstance != nullptr) {
169  AliError("An instance of AliAnalysisTaskEmcalEmbeddingHelper already exists: it will be deleted!!!");
170  delete fgInstance;
171  }
172 
173  fgInstance = this;
174 }
175 
182  AliAnalysisTaskSE(name),
183  fTriggerMask(0),
184  fMCRejectOutliers(false),
186  fZVertexCut(10),
187  fMaxVertexDist(999),
189  fInitializedNewFile(false),
190  fInitializedEmbedding(false),
191  fWrappedAroundTree(false),
192  fTreeName("aodTree"),
193  fNPtHardBins(1),
194  fPtHardBin(-1),
195  fRandomEventNumberAccess(kFALSE),
196  fRandomFileAccess(kTRUE),
197  fCreateHisto(true),
198  fYAMLConfig(),
202  fEmbeddedEventUsed(true),
203  fCentMin(-999),
204  fCentMax(-999),
206  fAutoConfigureBasePath("alien:///alice/cern.ch/user/a/alitrain/"),
207  fAutoConfigureTrainTypePath("PWGJE/Jets_EMC_PbPb/"),
208  fAutoConfigureIdentifier("autoConfigIdentifier"),
209  fFilePattern(""),
210  fInputFilename(""),
211  fFileListFilename(""),
212  fFilenameIndex(-1),
213  fFilenames(),
214  fConfigurationPath(""),
218  fChain(nullptr),
219  fCurrentEntry(0),
220  fLowerEntry(0),
221  fUpperEntry(0),
222  fOffset(0),
224  fFileNumber(0),
225  fHistManager(name),
226  fOutput(nullptr),
230  fPythiaTrials(0),
234  fPythiaPtHard(0.)
235 {
236  if (fgInstance != 0) {
237  AliError("An instance of AliAnalysisTaskEmcalEmbeddingHelper already exists: it will be deleted!!!");
238  delete fgInstance;
239  }
240 
241  fgInstance = this;
242 
243  if (fCreateHisto) {
244  DefineOutput(1, AliEmcalList::Class());
245  }
246 }
247 
254 {
255  if (fgInstance == this) fgInstance = nullptr;
256  if (fExternalEvent) delete fExternalEvent;
257  if (fExternalFile) {
258  fExternalFile->Close();
259  delete fExternalFile;
260  }
261 }
262 
264 {
265  // Initialize %YAML configuration, if one is given
266  bool initializedYAML = InitializeYamlConfig();
267 
269 
270  // Get file list
271  bool result = GetFilenames();
272 
273  if (result && initializedYAML) {
275  }
276 
277  if (removeDummyTask == true) {
278  RemoveDummyTask();
279  }
280 
281  // Initialize the YAML config object for streaming
283 
284  // Print the results of the initialization
285  // Print outside of the ALICE Log system to ensure that it is always available!
286  std::cout << *this;
287 
288  return result;
289 }
290 
298 {
299  // Following the variable blocks defined in the header
300  // Embedded event properties
301  std::vector<std::string> physicsSelection;
302  bool res = fYAMLConfig.GetProperty("embeddedEventPhysicsSelection", physicsSelection, false);
303  if (res) {
305  }
306  res = fYAMLConfig.GetProperty("enableMCOutlierRejection", fMCRejectOutliers, false);
307  res = fYAMLConfig.GetProperty("ptHardJetPtRejectionFactor", fPtHardJetPtRejectionFactor, false);
308  res = fYAMLConfig.GetProperty("embeddedEventZVertexCut", fZVertexCut, false);
309  res = fYAMLConfig.GetProperty("maxVertexDifferenceDistance", fMaxVertexDist, false);
310 
311  // Embedding helper properties
312  res = fYAMLConfig.GetProperty("treeName", fTreeName, false);
313  res = fYAMLConfig.GetProperty("nPtHardBins", fNPtHardBins, false);
314  res = fYAMLConfig.GetProperty("ptHardBin", fPtHardBin, false);
315  res = fYAMLConfig.GetProperty("randomEventNumberAccess", fRandomEventNumberAccess, false);
316  res = fYAMLConfig.GetProperty("randomFileAccess", fRandomFileAccess, false);
317  res = fYAMLConfig.GetProperty("createHisto", fCreateHisto, false);
318  // More general embedding helper properties
319  res = fYAMLConfig.GetProperty("filePattern", fFilePattern, false);
320  res = fYAMLConfig.GetProperty("inputFilename", fInputFilename, false);
321  res = fYAMLConfig.GetProperty("fileListFilename", fFileListFilename, false);
322  res = fYAMLConfig.GetProperty("filenameIndex", fFilenameIndex, false);
323  // Configuration path makes no sense, as we are already using the %YAML configuration
324  res = fYAMLConfig.GetProperty("runlist", fEmbeddedRunlist, false);
325  // Generally should not be set
326  res = fYAMLConfig.GetProperty("filenames", fFilenames, false);
327  res = fYAMLConfig.GetProperty("fPythiaCrossSectionFilenames", fPythiaCrossSectionFilenames, false);
328 
329  // Internal event selection properties
330  // NOTE: Need to define the base name here so that the property path is not ambiguous (due to otherwise only being `const char *`)
331  std::string baseName = "internalEventSelection";
332  res = fYAMLConfig.GetProperty({baseName, "enabled"}, fUseInternalEventSelection, false);
333  res = fYAMLConfig.GetProperty({baseName, "useManualCuts"}, fUseManualInternalEventCuts, false);
334  // Centrality
335  std::vector <double> centralityRange;
336  res = fYAMLConfig.GetProperty({baseName, "centralityRange"}, centralityRange, false);
337  if (res) {
338  if (centralityRange.size() != 2) {
339  AliErrorStream() << "Passed centrality range with " << centralityRange.size() << " entries, but 2 values are required. Ignoring values.\n";
340  }
341  else {
342  AliDebugStream(1) << "Setting internal event centrality range to [" << centralityRange.at(0) << ", " << centralityRange.at(1) << "]\n";
343  fCentMin = centralityRange.at(0);
344  fCentMax = centralityRange.at(1);
345  }
346  }
347  // Physics selection
348  res = fYAMLConfig.GetProperty({baseName, "physicsSelection"}, physicsSelection, false);
349  if (res) {
350  fOfflineTriggerMask = AliEmcalContainerUtils::DeterminePhysicsSelectionFromYAML(physicsSelection);
351  }
352 
353  // Auto configure pt hard properties
354  res = fYAMLConfig.GetProperty("autoConfigurePtHardBins", fAutoConfigurePtHardBins, false);
355  res = fYAMLConfig.GetProperty("autoConfigureBasePath", fAutoConfigureBasePath, false);
356  res = fYAMLConfig.GetProperty("autoConfigureTrainTypePath", fAutoConfigureTrainTypePath, false);
357  res = fYAMLConfig.GetProperty("autoConfigureIdentifier", fAutoConfigureIdentifier, false);
358 }
359 
384 {
385  // Determine the pattern filename if not yet set
386  if (fInputFilename == "") {
387  if (fTreeName == "aodTree") {
388  fInputFilename = "AliAOD.root";
389  }
390  else if (fTreeName == "esdTree") {
391  fInputFilename = "AliESDs.root";
392  }
393  else {
394  AliFatal(TString::Format("Requested default (pattern) input filename, but could not determine file type from tree name \"%s\". Please check the tree name and set the pattern input filename", fTreeName.Data()));
395  }
396  }
397 
398  // Retrieve filenames if we don't have them yet.
399  if (fFilenames.size() == 0)
400  {
401  // Handle pt hard bin auto configuration
403  {
404  if (fPtHardBin > 0) {
405  AliFatal("Requested both pt hard bin auto configuration and selected a non-zero pt hard bin. These are incompatible options. Please check your configuration.");
406  }
407  bool success = AutoConfigurePtHardBins();
408  if (success == false) {
409  AliFatal("Pt hard bin auto configuration requested, but it failed. Please check the logs.\n");
410  }
411  }
412 
413  // Handle if fPtHardBin is set
414  // This will require formatting the file pattern in the proper way to support these substitutions
415  if (fPtHardBin != -1 && fFilePattern != "") {
416  fFilePattern = TString::Format(fFilePattern, fPtHardBin);
417  }
418 
419  // Setup AliEn access if needed
420  if (fFilePattern.Contains("alien://") || fFileListFilename.Contains("alien://")) {
422  }
423 
424  // Retrieve AliEn filenames directly from AliEn
425  bool usedFilePattern = false;
426  if (fFilePattern.Contains("alien://")) {
427  usedFilePattern = true;
428  AliDebug(2, TString::Format("Trying to retrieve file list from AliEn with pattern file %s...", fFilePattern.Data()));
429 
430  // Create a temporary filename based on a UUID to make sure that it doesn't overwrite any files
431  if (fFileListFilename == "") {
433  }
434 
435  // The query command cannot handle "alien://" in the file pattern, so we need to remove it for the command
436  TString filePattern = fFilePattern;
437  filePattern.ReplaceAll("alien://", "");
438 
439  // Execute the grid query to get the filenames
440  AliDebug(2, TString::Format("Trying to retrieve file list from AliEn with pattern \"%s\" and input filename \"%s\"", filePattern.Data(), fInputFilename.Data()));
441  auto result = gGrid->Query(filePattern.Data(), fInputFilename.Data());
442 
443  if (result) {
444 
445  // Loop over the result to store it in the fileList file
446  std::ofstream outFile(fFileListFilename);
447  for (int i = 0; i < result->GetEntries(); i++)
448  {
449  TString path = result->GetKey(i, "turl");
450  // "turl" corresponds to the full AliEn url
451 
452  // If a runlist is specified for good embedded runs, only include the file if it is in this runlist
453  if (IsRunInRunlist(path.Data())) {
454  outFile << path << "\n";
455  }
456  }
457  outFile.close();
458  }
459  else {
460  AliErrorStream() << "Failed to run grid query\n";
461  return false;
462  }
463  }
464 
465  // Handle a filelist on AliEn
466  if (fFileListFilename.Contains("alien://")) {
467  // Check if we already used the file pattern
468  if (usedFilePattern) {
469  AliErrorStream() << "You set both the file pattern and the file list filename! The file list filename will override the pattern! Pattern: \"" << fFilePattern << "\", filename: \"" << fFileListFilename << "\"\nPlease check that this is the desired behavior!\n";
470  }
471 
472  // Determine the local filename and copy file to local directory
473  std::string alienFilename = fFileListFilename.Data();
474  fFileListFilename = gSystem->BaseName(alienFilename.c_str());
476 
477  TFile::Cp(alienFilename.c_str(), fFileListFilename.Data());
478  }
479 
480  std::ifstream inputFile(fFileListFilename);
481 
482  // Copy available files into the filenames vector
483  // From:: https://stackoverflow.com/a/8365247
484  std::copy(std::istream_iterator<std::string>(inputFile),
485  std::istream_iterator<std::string>(),
486  std::back_inserter(fFilenames));
487 
488  inputFile.close();
489  }
490 
491  if (fFilenames.size() == 0) {
492  AliFatal(TString::Format("Filenames from pattern \"%s\" and file list \"%s\" yielded an empty list!", fFilePattern.Data(), fFileListFilename.Data()));
493  }
494 
495  // Add "#" to files in there are any zip files
496  // It is require to open the proper root file within the zip
497  for (auto filename : fFilenames)
498  {
499  if (filename.find(".zip") != std::string::npos && filename.find("#") == std::string::npos) {
500  filename += "#";
502  if (fTreeName == "aodTree") {
503  filename += "#AliAOD.root";
504  }
505  else if (fTreeName == "esdTree") {
506  filename += "#AliESDs.root";
507  }
508  else {
509  AliError(TString::Format("Filename %s contains \".zip\" and not \"#\", but tree name %s is not recognized. Please check the file list to ensure that the proper path is set.", filename.c_str(), fTreeName.Data()));
510  return false;
511  }
512  }
513  }
514 
515  // Determine whether AliEn is needed
516  // It is possible that this has not been determined up to this point
517  for (auto filename : fFilenames)
518  {
519  if (filename.find("alien://") != std::string::npos) {
521  // No point in continuing to search once we know that it is needed
522  break;
523  }
524  }
525 
526  // Check if each filenames exists. If they do not exist, then remove them for fFilenames
527  unsigned int initialSize = fFilenames.size();
528  // NOTE: We invert the result of IsFileAccessible because we should return true for files that should be _removed_ (ie are inaccessible)
529  fFilenames.erase(std::remove_if(fFilenames.begin(), fFilenames.end(), [](const std::string & filename) {return (::IsFileAccessible(filename) == false);} ), fFilenames.end());
530 
531  AliInfoStream() << "Found " << fFilenames.size() << " files to embed (" << (initialSize - fFilenames.size()) << " filename(s) inaccessible or invalid)\n";
532 
533  // Determine pythia filename
535 
536  return true;
537 }
538 
544 {
545  // Get the initial filename. Use the first entry as a proxy for other input files
546  std::string externalEventFilename = "";
547  if (fFilenames.size() > 0) {
548  externalEventFilename = fFilenames.at(0);
549  }
550  else {
551  return;
552  }
553 
554  std::vector <std::string> pythiaBaseFilenames = {"pyxsec.root", "pyxsec_hists.root"};
555  AliInfoStream() << "Attempting to determine pythia cross section filename. It can be normal to see some TFile::Init() errors!\n";
556  std::string pythiaXSecFilename = "";
557  for (auto & name : pythiaBaseFilenames) {
558  pythiaXSecFilename = ConstructFullPythiaXSecFilename(externalEventFilename, name, true);
559  if (pythiaXSecFilename != "") {
560  AliDebugStream(4) << "Found pythia cross section filename \"" << name << "\"\n";
561  fPythiaXSecFilename = name;
562  break;
563  }
564  }
565 
566  if (fPythiaXSecFilename == "") {
567  // Failed entirely - just give up on this
568  // We will use an empty filename as a proxy for whether the file has been found (empty is equivalent to not found)
569  AliErrorStream() << "Failed to find pythia x sec file! Continuing with only the pythia header!\n";
570  }
571  else {
572  AliInfoStream() << "Found pythia cross section file \"" << fPythiaXSecFilename << "\".\n";
573  }
574 }
575 
576 
584 bool AliAnalysisTaskEmcalEmbeddingHelper::IsRunInRunlist(const std::string & path) const
585 {
586  if (fEmbeddedRunlist.size() == 0) {
587  return true;
588  }
589 
590  for (auto run : fEmbeddedRunlist) {
591  if (path.find(run) != std::string::npos) {
592  return true;
593  }
594  }
595  return false;
596 }
597 
604 {
605  if (fConfigurationPath == "") {
606  AliInfo("No Embedding YAML configuration was provided");
607  }
608  else {
609  AliInfoStream() << "Embedding YAML configuration was provided: \"" << fConfigurationPath << "\".\n";
610 
611  int addedConfig = fYAMLConfig.AddConfiguration(fConfigurationPath, "yamlConfig");
612  if (addedConfig < 0) {
613  AliError("YAML Configuration not found!");
614  return false;
615  }
616  }
617 
618  return true;
619 }
620 
637 {
638  bool returnValue = false;
639 
640  AliInfoStream() << "Attempting to auto configure pt hard bins.\n";
641  // %YAML configuration containing pt hard bin to train number mapping
643 
644  // Handle AliEn explicitly here since the default base path contains "alien://"
645  if (fAutoConfigureBasePath.find("alien://") != std::string::npos && !gGrid) {
647  }
648 
649  // Get train ID
650  // Need to get the char * directly because it may be null.
651  const char * trainNumberStr = gSystem->Getenv("TRAIN_RUN_ID");
652  std::stringstream trainNumberSS;
653  if (trainNumberStr) {
654  trainNumberSS << trainNumberStr;
655  }
656  if (trainNumberSS.str() == "") {
657  AliFatal("Cannot retrieve train ID.");
658  }
659  // Extract train number from the string
660  int trainNumber;
661  trainNumberSS >> trainNumber;
662 
663  // Determine the file path
665  filename += "/";
667  filename += "/";
669  // Add ".yaml" if it is not already there
670  std::string yamlExtension = ".yaml";
671  if (filename.find(yamlExtension) == std::string::npos) {
672  filename += yamlExtension;
673  }
674 
675  // Check if file exists
676  if (gSystem->AccessPathName(filename.c_str())) {
677  // File _does not_ exist
678  AliInfoStream() << "Train pt hard bin configuration file not available, so creating a new empty configuration named \"" << fAutoConfigureIdentifier << "\".\n";
679  // Use an empty configuration
681  }
682  else {
683  AliInfoStream() << "Opening configuration located at \"" << filename << "\".\n";
684  // Use the existing configuration
686  }
687 
688  // Look for each pt hard bin, and then retrieve the corresponding train number
689  // Once an open pt hard bin is found, add the current train number
690  int tempTrainNumber = -1;
691  bool getPropertyReturnValue = false;
692  std::stringstream propertyName;
693  for (int ptHardBin = 1; ptHardBin <= fNPtHardBins; ptHardBin++)
694  {
695  propertyName.str("");
696  propertyName << ptHardBin;
697  getPropertyReturnValue = config.GetProperty(propertyName.str(), tempTrainNumber, false);
698  if (getPropertyReturnValue != true) {
699  AliInfoStream() << "Train " << trainNumber << " will use pt hard bin " << ptHardBin << ".\n";
700  // We have determine our pt hard bin!
701  fPtHardBin = ptHardBin;
702 
703  // Write the train number back out to the %YAML configuration and save it
704  config.WriteProperty(propertyName.str(), trainNumber, fAutoConfigureIdentifier);
706 
707  // NOTE: Cannot clean up the yaml file on the last pt hard bin because the train can be launched
708  // multiple times due to tests, etc. Therefore, we have to accept that we are leaving around used
709  // yaml config files.
710 
711  // We are done - continue on.
712  returnValue = true;
713  break;
714  }
715  else {
716  AliDebugStream(2) << "Found pt hard bin " << ptHardBin << " corresponding to train number " << trainNumber << ".\n";
717  // If train was already allocated (say, by a test train), then use that pt hard bin
718  if (tempTrainNumber == trainNumber) {
719  AliInfoStream() << "Train run number " << trainNumber << " was already found assigned to pt hard bin " << ptHardBin << ". That pt hard bin will be used.\n";
720  fPtHardBin = ptHardBin;
721 
722  // We are done - continue on.
723  returnValue = true;
724  break;
725  }
726  // Otherwise, nothing to be done.
727  }
728  }
729 
730  return returnValue;
731 }
732 
741 {
742  std::string tempStr = "";
743  if (fFileListFilename == "") {
744  tempStr = "fileList";
745  }
746  TUUID tempUUID;
747  tempStr += ".";
748  tempStr += tempUUID.AsString();
749  tempStr += ".txt";
750 
751  return tempStr;
752 }
753 
762 {
763  while (filename.rbegin() != filename.rend() && *(filename.rbegin()) == '/') {
764  filename.pop_back();
765  }
766 
767  return filename;
768 }
769 
775 {
776  // This determines which file is added first to the TChain, thus determining the order of processing
777  // Random file access. Only do this if the user has no set the filename index and request random file access
778  if (fFilenameIndex == -1 && fRandomFileAccess) {
779  // Floor ensures that we it doesn't overflow
780  TRandom3 rand(0);
781  fFilenameIndex = TMath::FloorNint(rand.Rndm()*fFilenames.size());
782  // +1 to account for the fact that the filenames vector is 0 indexed.
783  AliInfo(TString::Format("Starting with random file number %i!", fFilenameIndex+1));
784  }
785  // If not random file access, then start from the beginning
786  if (fFilenameIndex < 0 || static_cast<UInt_t>(fFilenameIndex) >= fFilenames.size()) {
787  // Skip notifying on -1 since it will likely be set there due to constructor.
788  if (fFilenameIndex != -1) {
789  AliWarning(TString::Format("File index %i out of range from 0 to %lu! Resetting to 0!", fFilenameIndex, fFilenames.size()));
790  }
791  fFilenameIndex = 0;
792  }
793 
794  // +1 to account for the fact that the filenames vector is 0 indexed.
795  AliInfo(TString::Format("Starting with file number %i out of %lu", fFilenameIndex+1, fFilenames.size()));
796 }
797 
807 {
808  Int_t attempts = -1;
809 
810  do {
811  // Reset to start of tree
812  if (fCurrentEntry == fUpperEntry) {
814  fWrappedAroundTree = true;
815  }
816 
818  // Continue with GetEntry as normal
819  }
820  else {
821  // NOTE: On transition from one file to the next, this calls the next entry that would be expected.
822  // However, if it is for the last file, it tries to GetEntry() of one entry past the end of the last file.
823  // Normally, this would be a problem, however GetEntry() just doesn't fill the fields of an invalid index
824  // instead of throwing an error. So "invalid values" are filled for a file that doesn't exist, but then
825  // they are immediately replaced by the lines below that reset the access values and re-init the tree.
826  // The benefit of this approach is it simplies file counting (we don't need to carefully increment here
827  // and in InitTree()) and preserves the desired behavior when we are not at the last file.
828  InitTree();
829  }
830 
831  // Load current event
832  // Can be a simple less than, because fFileNumber counts from 0.
834  fChain->GetEntry(fCurrentEntry);
835  }
836  else {
837  AliError("====================================================================================================");
838  AliError("== No more files available to embed from the TChain! Restarting from the beginning of the TChain! ==");
839  AliError("== Be careful to check that this is the desired action! ==");
840  AliError("====================================================================================================");
841 
842  // Reset the relevant access values
843  // fCurrentEntry and fLowerEntry are automatically reset in InitTree()
844  fFileNumber = 0;
845  fUpperEntry = 0;
846 
847  // Re-init back to the start
848  InitTree();
849 
850  // Access the relevant entry
851  // We are certain that fFileNumber is less than fMaxNumberOfFiles, so we are resetting to start
852  fChain->GetEntry(fCurrentEntry);
853  }
854  AliDebug(4, TString::Format("Loading entry %i between %i-%i, starting with offset %i from the lower bound of %i", fCurrentEntry, fLowerEntry, fUpperEntry, fOffset, fLowerEntry));
855 
856  // Set relevant event properties
858 
859  // Increment current entry
860  fCurrentEntry++;
861 
862  // Provide a check for number of attempts
863  attempts++;
864  if (attempts == 1000)
865  AliWarning("After 1000 attempts no event has been accepted by the event selection (trigger, centrality...)!");
866 
867  // Record event properties
868  if (fCreateHisto) {
870  }
871 
872  } while (!IsEventSelected());
873 
874  if (fCreateHisto) {
875  fHistManager.FillTH1("fHistEventCount", "Accepted");
876  fHistManager.FillTH1("fHistEmbeddedEventsAttempted", attempts);
877  }
878 
879  if (!fChain) return kFALSE;
880 
881  return kTRUE;
882 }
883 
889 {
890  AliDebug(4, "Set event properties");
891  fExternalHeader = fExternalEvent->GetHeader();
892 
893  // Handle pythia header if AOD
894  AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(fExternalEvent->FindListObject(AliAODMCHeader::StdBranchName()));
895  if (aodMCH) {
896  for (UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++) {
897  fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
898  if (fPythiaHeader) break;
899  }
900  }
901 
902  if (fPythiaHeader)
903  {
904  fPythiaCrossSection = fPythiaHeader->GetXsection();
905  fPythiaTrials = fPythiaHeader->Trials();
906  fPythiaPtHard = fPythiaHeader->GetPtHard();
907  // It is identically zero if the cross section is not available
908  if (fPythiaCrossSection == 0.) {
909  AliDebugStream(4) << "Taking the pythia cross section avg from the xsec file.\n";
911  }
912  // It is identically zero if the number of trials is not available
913  if (fPythiaTrials == 0.) {
914  AliDebugStream(4) << "Taking the pythia trials avg from the xsec file.\n";
916  }
917  // Pt hard is inherently event-by-event and cannot by taken as a avg quantity.
918 
919  AliDebugStream(4) << "Pythia header is defined!\n";
920  AliDebugStream(4) << "fPythiaCrossSection: " << fPythiaCrossSection << "\n";
921  }
922 }
923 
928 {
929  // Fill trials, xsec, pt hard
930  fHistManager.FillTH1("fHistTrials", fPtHardBin, fPythiaTrials);
932  fHistManager.FillTH1("fHistPtHard", fPythiaPtHard);
933 }
934 
941 {
943  return kTRUE;
944  }
945 
946  if (fCreateHisto) {
947  // Keep count of number of rejected events
948  fHistManager.FillTH1("fHistEventCount", "Rejected");
949  }
950 
951  return kFALSE;
952 }
953 
960 {
961  // Check if pt hard bin is 0, indicating a problem with the event or the grid.
962  // In such a case, the event should be rejected.
963  // This condition should only be applied if we have a valid pythia header.
964  // (pt hard should still be set even if the production wasn't done in pt hard bins).
965  if (fPythiaPtHard == 0. && fPythiaHeader) {
966  AliDebugStream(3) << "Event rejected due to pt hard = 0, indicating a problem with the external event.\n";
967  if (fCreateHisto) {
968  fHistManager.FillTH1("fHistEmbeddedEventRejection", "PtHardIs0", 1);
969  }
970  return kFALSE;
971  }
972 
973  // Physics selection
974  if (fTriggerMask != 0) {
975  UInt_t res = 0;
976  const AliESDEvent *eev = dynamic_cast<const AliESDEvent*>(fExternalEvent);
977  if (eev) {
978  AliFatal("Event selection is not implemented for embedding ESDs.");
979  // Unfortunately, the normal method of retrieving the trigger mask (commented out below) doesn't work for the embedded event since we don't
980  // create an input handler and I am not an expert on getting a trigger mask. Further, embedding ESDs is likely to be inefficient, so it is
981  // probably best to avoid it if possible.
982  //
983  // Suggestions are welcome here!
984  //res = (dynamic_cast<AliInputEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
985  } else {
986  const AliAODEvent *aev = dynamic_cast<const AliAODEvent*>(fExternalEvent);
987  if (aev) {
988  res = (dynamic_cast<AliVAODHeader*>(aev->GetHeader()))->GetOfflineTrigger();
989  }
990  }
991 
992  if ((res & fTriggerMask) == 0) {
993  AliDebug(3, Form("Event rejected due to physics selection. Event trigger mask: %d, trigger mask selection: %d.",
994  res, fTriggerMask));
995  if (fCreateHisto) {
996  fHistManager.FillTH1("fHistEmbeddedEventRejection", "PhysSel", 1);
997  }
998  return kFALSE;
999  }
1000  }
1001 
1002  // Vertex selection
1003  Double_t externalVertex[3]={0};
1004  Double_t inputVertex[3]={0};
1005  const AliVVertex *externalVert = fExternalEvent->GetPrimaryVertex();
1006  const AliVVertex *inputVert = AliAnalysisTaskSE::InputEvent()->GetPrimaryVertex();
1007  if (externalVert && inputVert) {
1008  externalVert->GetXYZ(externalVertex);
1009  inputVert->GetXYZ(inputVertex);
1010 
1011  if (TMath::Abs(externalVertex[2]) > fZVertexCut) {
1012  AliDebug(3, Form("Event rejected due to Z vertex selection. Event Z vertex: %f, Z vertex cut: %f",
1013  externalVertex[2], fZVertexCut));
1014  if (fCreateHisto) {
1015  fHistManager.FillTH1("fHistEmbeddedEventRejection", "Vz", 1);
1016  }
1017  return kFALSE;
1018  }
1019  Double_t dist = TMath::Sqrt((externalVertex[0]-inputVertex[0])*(externalVertex[0]-inputVertex[0])+(externalVertex[1]-inputVertex[1])*(externalVertex[1]-inputVertex[1])+(externalVertex[2]-inputVertex[2])*(externalVertex[2]-inputVertex[2]));
1020  if (dist > fMaxVertexDist) {
1021  AliDebug(3, Form("Event rejected because the distance between the current and embedded vertices is > %f. "
1022  "Current event vertex (%f, %f, %f), embedded event vertex (%f, %f, %f). Distance = %f",
1023  fMaxVertexDist, inputVertex[0], inputVertex[1], inputVertex[2], externalVertex[0], externalVertex[1], externalVertex[2], dist));
1024  if (fCreateHisto) {
1025  fHistManager.FillTH1("fHistEmbeddedEventRejection", "VertexDist", 1);
1026  }
1027  return kFALSE;
1028  }
1029  }
1030 
1031  // Check for pt hard bin outliers
1033  {
1034  // Pythia jet / pT-hard > factor
1035  // This corresponds to "condition 1" in AliAnalysisTaskEmcal
1036  // NOTE: The other "conditions" defined there are not really suitable to define here, since they
1037  // depend on the input objects of the event
1038  if (fPtHardJetPtRejectionFactor > 0.) {
1039  TLorentzVector jet;
1040 
1041  Int_t nTriggerJets = fPythiaHeader->NTriggerJets();
1042 
1043  AliDebugStream(4) << "Pythia Njets: " << nTriggerJets << ", pT Hard: " << fPythiaPtHard << "\n";
1044 
1045  Float_t tmpjet[]={0,0,0,0};
1046  for (Int_t iJet = 0; iJet< nTriggerJets; iJet++) {
1047  fPythiaHeader->TriggerJet(iJet, tmpjet);
1048 
1049  jet.SetPxPyPzE(tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3]);
1050 
1051  AliDebugStream(5) << "Pythia jet " << iJet << ", pycell jet pT: " << jet.Pt() << "\n";
1052 
1053  //Compare jet pT and pt Hard
1054  if (jet.Pt() > fPtHardJetPtRejectionFactor * fPythiaPtHard) {
1055  AliDebugStream(3) << "Event rejected because of MC outlier removal. Pythia header jet with: pT Hard " << fPythiaPtHard << ", pycell jet pT " << jet.Pt() << ", rejection factor " << fPtHardJetPtRejectionFactor << "\n";
1056  fHistManager.FillTH1("fHistEmbeddedEventRejection", "MCOutlier", 1);
1057  return kFALSE;
1058  }
1059  }
1060  }
1061  }
1062 
1063  return kTRUE;
1064 }
1065 
1072 {
1073  if (!fChain) return kFALSE;
1074 
1075  if (!fExternalEvent) {
1076  if (fTreeName == "aodTree") {
1077  fExternalEvent = new AliAODEvent();
1078  }
1079  else if (fTreeName == "esdTree") {
1080  fExternalEvent = new AliESDEvent();
1081  }
1082  else {
1083  AliError(Form("Tree name %s not recognized!", fTreeName.Data()));
1084  return kFALSE;
1085  }
1086  }
1087 
1088  fExternalEvent->ReadFromTree(fChain, fTreeName);
1089 
1090  return kTRUE;
1091 }
1092 
1099 {
1100  SetupEmbedding();
1101 
1102  // Reinitialize the YAML config after it was streamed so that it can be used properly.
1104 
1105  // Setup AliEventCuts
1107  {
1108  AliDebugStream(1) << "Configuring AliEventCuts for internal event selection.\n";
1109  // Handle manual cuts
1111  fInternalEventCuts.SetManualMode();
1112  // Implement these cuts by retrieving the event cuts object and setting them manually.
1113  }
1114 
1115  // Trigger selection
1116  bool useEventCutsAutomaticTriggerSelection = false;
1117  bool res = fYAMLConfig.GetProperty(std::vector<std::string>({"internalEventSelection", "useEventCutsAutomaticTriggerSelection"}), useEventCutsAutomaticTriggerSelection, false);
1118  if (res && useEventCutsAutomaticTriggerSelection) {
1119  // Use the autmoatic selection. Nothing to be done.
1120  AliDebugStream(1) << "Using the automatic trigger selection from AliEventCuts.\n";
1121  }
1122  else {
1123  // Use the cuts selected by SelectCollisionCandidates()
1124  AliDebugStream(1) << "Using the trigger selection specified with SelectCollisionCandidates()\n.";
1125  fInternalEventCuts.OverrideAutomaticTriggerSelection(fOfflineTriggerMask);
1126  }
1127  }
1128 
1129  if (!fCreateHisto) {
1130  return;
1131  }
1132 
1133  // Create output list
1134  OpenFile(1);
1135  fOutput = new AliEmcalList();
1136  fOutput->SetOwner();
1137 
1138  // Get the histograms from AliEventCuts
1140  // This list will be owned by fOutput, so it won't be leaked.
1141  TList * eventCutsOutput = new TList();
1142  eventCutsOutput->SetOwner(kTRUE);
1143  eventCutsOutput->SetName("EventCuts");
1144 
1145  // Add the event cuts to the output
1146  fInternalEventCuts.AddQAplotsToList(eventCutsOutput);
1147  fOutput->Add(eventCutsOutput);
1148  }
1149 
1150  // Create histograms
1151  TString histName;
1152  TString histTitle;
1153 
1154  // Cross section
1155  histName = "fHistXsection";
1156  histTitle = "Pythia Cross Section;p_{T} hard bin; XSection";
1157  fHistManager.CreateTProfile(histName, histTitle, fNPtHardBins, 0, fNPtHardBins);
1158 
1159  // Trials
1160  histName = "fHistTrials";
1161  histTitle = "Number of Pythia Trials;p_{T} hard bin;Trials";
1162  fHistManager.CreateTH1(histName, histTitle, fNPtHardBins, 0, fNPtHardBins);
1163 
1164  // Pt hard spectra
1165  histName = "fHistPtHard";
1166  histTitle = "p_{T} Hard Spectra;p_{T} hard;Counts";
1167  fHistManager.CreateTH1(histName, histTitle, 500, 0, 1000);
1168 
1169  // Count of accepted and rejected events
1170  histName = "fHistEventCount";
1171  histTitle = "Event count;Result;Count";
1172  auto histEventCount = fHistManager.CreateTH1(histName, histTitle, 2, 0, 2);
1173  histEventCount->GetXaxis()->SetBinLabel(1,"Accepted");
1174  histEventCount->GetXaxis()->SetBinLabel(2,"Rejected");
1175 
1176  // Event rejection reason
1177  histName = "fHistEmbeddedEventRejection";
1178  histTitle = "Reasons to reject embedded event";
1179  std::vector<std::string> binLabels = {"PhysSel", "MCOutlier", "Vz", "VertexDist", "PtHardIs0"};
1180  auto histEmbeddedEventRejection = fHistManager.CreateTH1(histName, histTitle, binLabels.size(), 0, binLabels.size());
1181  // Set label names
1182  for (unsigned int i = 1; i <= binLabels.size(); i++) {
1183  histEmbeddedEventRejection->GetXaxis()->SetBinLabel(i, binLabels.at(i-1).c_str());
1184  }
1185  histEmbeddedEventRejection->GetYaxis()->SetTitle("Counts");
1186 
1187  // Rejected events in embedded event selection
1188  histName = "fHistEmbeddedEventsAttempted";
1189  histTitle = "Number of embedded events rejected by event selection before success;Number of rejected events;Counts";
1190  fHistManager.CreateTH1(histName, histTitle, 200, 0, 200);
1191 
1192  // Number of files embedded
1193  histName = "fHistNumberOfFilesEmbedded";
1194  histTitle = "Number of files which contributed events to be embedded";
1195  fHistManager.CreateTH1(histName, histTitle, 1, 0, 2);
1196 
1197  // File number which was embedded
1198  histName = "fHistAbsoluteFileNumber";
1199  histTitle = "Number of times each absolute file number was embedded";
1200  fHistManager.CreateTH1(histName, histTitle, fMaxNumberOfFiles, 0, fMaxNumberOfFiles);
1201 
1203  // Internal event cut statistics
1204  histName = "fHistInternalEventCutsStats";
1205  histTitle = "Number of events to pass each cut";
1206  binLabels = {"passedEventCuts", "centrality", "passedAllCuts"};
1207  auto histInternalEventCutsStats = fHistManager.CreateTH1(histName, histTitle, binLabels.size(), 0, binLabels.size());
1208  // Set label names
1209  for (unsigned int i = 1; i <= binLabels.size(); i++) {
1210  histInternalEventCutsStats->GetXaxis()->SetBinLabel(i, binLabels.at(i-1).c_str());
1211  }
1212  histInternalEventCutsStats->GetYaxis()->SetTitle("Number of selected events");
1213  }
1214 
1215  // Add all histograms to output list
1216  TIter next(fHistManager.GetListOfHistograms());
1217  TObject* obj = 0;
1218  while ((obj = next())) {
1219  fOutput->Add(obj);
1220  }
1221 
1222  PostData(1, fOutput);
1223 }
1224 
1232 {
1233  // Determine which file to start with
1235 
1236  // Setup TChain
1237  fChain = new TChain(fTreeName);
1238 
1239  // Determine whether AliEn is needed
1240  for (auto filename : fFilenames)
1241  {
1242  if (filename.find("alien://") != std::string::npos) {
1243  ::ConnectToAliEn();
1244  // No point in continuing to search once we know that it is needed
1245  break;
1246  }
1247  }
1248 
1249  // Add files for TChain
1250  // See: https://stackoverflow.com/a/8533198
1251  bool wrapped = false;
1252  std::string fullPythiaXSecFilename = "";
1253  for (auto filename = fFilenames.begin() + fFilenameIndex; (filename != fFilenames.begin() + fFilenameIndex || !wrapped); filename++)
1254  {
1255  // Wraps the loop back around to the beginning
1256  if (filename == fFilenames.end())
1257  {
1258  // Explicit check is needed. Otherwise, an offset of 0 would load the 0th entry twice.
1259  if (fFilenameIndex == 0) {
1260  break;
1261  }
1262  filename = fFilenames.begin();
1263  wrapped = true;
1264  }
1265 
1266  // Add to the Chain
1267  AliDebugStream(4) << "Adding file to the embedded input chain \"" << *filename << "\".\n";
1268  fChain->Add(filename->c_str());
1269 
1270  // Determine the full pythia cross section filename based on the previously determined filename
1271  // If we have determined that it doesn't exist in the initialization then we don't repeated attempt to open
1272  // the file (which will fail)
1273  if (fPythiaXSecFilename != "") {
1274  // Could check here again whether it exists here, but almost certainly unnecessary.
1275  // Further, we won't check to ensure that rapid, repeated file access on AliEn doesn't cause any problmes!
1276  fullPythiaXSecFilename = ConstructFullPythiaXSecFilename(*filename, fPythiaXSecFilename, false);
1277 
1278  AliDebugStream(4) << "Adding pythia cross section file \"" << fullPythiaXSecFilename << "\".\n";
1279 
1280  // They will automatically be ordered the same as the files to embed!
1281  fPythiaCrossSectionFilenames.push_back(fullPythiaXSecFilename);
1282  }
1283  }
1284 
1285  // Keep track of the total number of files in the TChain to ensure that we don't start repeating within the chain
1286  fMaxNumberOfFiles = fChain->GetListOfFiles()->GetEntries();
1287 
1288  if (fFilenames.size() > fMaxNumberOfFiles) {
1289  AliErrorStream() << "Number of input files (" << fFilenames.size() << ") is larger than the number of available files (" << fMaxNumberOfFiles << "). Something went wrong when adding some of those files to the TChain!\n";
1290  }
1291 
1292  // Setup input event
1293  Bool_t res = InitEvent();
1294  if (!res) return kFALSE;
1295 
1296  return kTRUE;
1297 }
1298 
1309 std::string AliAnalysisTaskEmcalEmbeddingHelper::ConstructFullPythiaXSecFilename(std::string externalEventFilename, const std::string & pythiaFilename, bool testIfExists) const
1310 {
1311  std::string pythiaXSecFilename = "";
1312 
1313  // Remove "#*.root" if necessary
1314  if (externalEventFilename.find(".zip#") != std::string::npos) {
1315  std::size_t pos = externalEventFilename.find_last_of("#");
1316  externalEventFilename.erase(pos);
1317  }
1318 
1319  // Handle different file types
1320  if (externalEventFilename.find(".zip") != std::string::npos)
1321  {
1322  // Handle zip files
1323  pythiaXSecFilename = externalEventFilename;
1324  pythiaXSecFilename += "#";
1325  pythiaXSecFilename += pythiaFilename;
1326 
1327  // Check if the file is accessible
1328  if (testIfExists) {
1329  // Unfortunately, we cannot test for the existence of a file in an archive.
1330  // Instead, we have to tolerate TFile throwing an error (maximum of two).
1331  std::unique_ptr<TFile> fTemp(TFile::Open(pythiaXSecFilename.c_str(), "READ"));
1332 
1333  if (!fTemp) {
1334  AliDebugStream(4) << "File " << pythiaXSecFilename.c_str() << " does not exist!\n";
1335  pythiaXSecFilename = "";
1336  }
1337  else {
1338  AliDebugStream(4) << "Found pythia cross section file \"" << pythiaXSecFilename.c_str() << "\".\n";
1339  }
1340  }
1341  }
1342  else
1343  {
1344  // Handle normal root files
1345  pythiaXSecFilename = gSystem->DirName(externalEventFilename.c_str());
1346  pythiaXSecFilename += "/";
1347  pythiaXSecFilename += pythiaFilename;
1348 
1349  // Check if the file is accessible
1350  if (testIfExists) {
1351  if(::IsFileAccessible(pythiaXSecFilename)) {
1352  AliDebugStream(4) << "Found pythia cross section file \"" << pythiaXSecFilename.c_str() << "\".\n";
1353  }
1354  else {
1355  AliDebugStream(4) << "File " << pythiaXSecFilename.c_str() << " does not exist!\n";
1356  pythiaXSecFilename = "";
1357  }
1358  }
1359  }
1360 
1361  return pythiaXSecFilename;
1362 }
1363 
1370 {
1371  if (fInitializedConfiguration == false) {
1372  AliFatal("The configuration is not initialized. Check that Initialize() was called!");
1373  }
1374 
1375  // Setup TChain
1376  Bool_t res = SetupInputFiles();
1377  if (!res) { return; }
1378 
1379  // Note if getting random event access
1381  AliInfo("Random event number access enabled!");
1382  }
1383 
1384  fInitializedEmbedding = kTRUE;
1385 }
1386 
1398 {
1399  // Load first entry of the (next) file so that we can query information about it
1400  // (it is unaccessible otherwise).
1401  // Since fUpperEntry is the total number of entries, loading it will retrieve the
1402  // next tree (in the next file) since entries are indexed starting from 0.
1403  fChain->GetEntry(fUpperEntry);
1404 
1405  // Determine tree size and current entry
1406  // Set the limits of the new tree
1408  // Fine to be += as long as we started at 0
1409  fUpperEntry += fChain->GetTree()->GetEntries();
1410 
1411  // Jump ahead at random if desired
1412  // Determines the offset into the tree
1414  TRandom3 rand(0);
1415  fOffset = TMath::Nint(rand.Rndm()*(fUpperEntry-fLowerEntry))-1;
1416  }
1417  else {
1418  fOffset = 0;
1419  }
1420 
1421  // Sets which entry to start if the try
1423 
1424  // Keep track of the number of files that we have gone through
1425  // To start from 0, we only increment if fLowerEntry > 0
1426  if (fLowerEntry > 0) {
1427  fFileNumber++;
1428  }
1429 
1430  // Add to the count the number of files which were embedded
1431  fHistManager.FillTH1("fHistNumberOfFilesEmbedded", 1);
1432  fHistManager.FillTH1("fHistAbsoluteFileNumber", (fFileNumber + fFilenameIndex) % fMaxNumberOfFiles);
1433 
1434  // Check for pythia cross section and extract if possible
1435  // fFileNumber corresponds to the next file
1436  // If there are pythia filenames, the number of match the file number of the tree.
1437  // If we previously gave up on extracting then there should be no entires
1438  if (fPythiaCrossSectionFilenames.size() > 0) {
1439  // Need to check that fFileNumber is smaller than the size of the vector because we don't check if
1442 
1443  if (!success) {
1444  AliDebugStream(3) << "Failed to retrieve cross section from xsec file. Will still attempt to get the information from the header.\n";
1445  }
1446  }
1447  else {
1448  AliErrorStream() << "Attempted to read past the end of the pythia cross section filenames vector. File number: " << fFileNumber << ", vector size: " << fPythiaCrossSectionFilenames.size() << ".\nThis should only occur if we have run out of files to embed!\n";
1449  }
1450  }
1451 
1452  AliDebug(2, TString::Format("Will start embedding file %i beginning from entry %i (entry %i within the file). NOTE: This file number is not equal to the absolute file number in the file list!", fFileNumber, fCurrentEntry, fCurrentEntry - fLowerEntry));
1453  // NOTE: Cannot use this print message, as it is possible that fMaxNumberOfFiles != fFilenames.size() because
1454  // invalid filenames may be included in the fFilenames count!
1455  //AliDebug(2, TString::Format("Will start embedding file %i as the %ith file beginning from entry %i.", (fFilenameIndex + fFileNumber) % fMaxNumberOfFiles, fFileNumber, fCurrentEntry));
1456 
1457  // (re)set whether we have wrapped the tree
1458  fWrappedAroundTree = false;
1459 
1460  // Note that the tree in the new file has been initialized
1461  fInitializedNewFile = kTRUE;
1462 }
1463 
1472 {
1473  std::unique_ptr<TFile> fxsec(TFile::Open(pythiaFileName.c_str()));
1474 
1475  if (fxsec)
1476  {
1477  int trials = 0;
1478  double crossSection = 0;
1479  double nEvents = 0;
1480  // Check if it's a tree
1481  TTree *xtree = dynamic_cast<TTree*>(fxsec->Get("Xsection"));
1482  if (xtree) {
1483  UInt_t ntrials = 0;
1484  Double_t xsection = 0;
1485  xtree->SetBranchAddress("xsection",&xsection);
1486  xtree->SetBranchAddress("ntrials",&ntrials);
1487  xtree->GetEntry(0);
1488  trials = ntrials;
1489  crossSection = xsection;
1490  // TODO: Test this on a file which has pyxsec.root!
1491  nEvents = 1.;
1492  AliFatal("Have no tested pyxsec.root files. Need to determine the proper way to get nevents!!");
1493  }
1494  else {
1495  // Check if it's instead the histograms
1496  // find the tlist we want to be independtent of the name so use the Tkey
1497  TKey* key = static_cast<TKey*>(fxsec->GetListOfKeys()->At(0));
1498  if (!key) return false;
1499  TList *list = dynamic_cast<TList*>(key->ReadObj());
1500  if (!list) return false;
1501  TProfile * crossSectionHist = static_cast<TProfile*>(list->FindObject("h1Xsec"));
1502  // check for failure
1503  if(!(crossSectionHist->GetEntries())) {
1504  // No cross seciton information available - fall back to raw
1505  AliErrorStream() << "No cross section information available in file \"" << fxsec->GetName() << "\". Will still attempt to extract cross section information from pythia header.\n";
1506  } else {
1507  // Cross section histogram filled - take it from there
1508  crossSection = crossSectionHist->GetBinContent(1);
1509  if(!crossSection) AliErrorStream() << GetName() << ": Cross section 0 for file " << pythiaFileName << std::endl;
1510  }
1511  TH1F * trialsHist = static_cast<TH1F*>(list->FindObject("h1Trials"));
1512  trials = trialsHist->GetBinContent(1);
1513  nEvents = trialsHist->GetEntries();
1514  }
1515 
1516  // If successful in retrieveing the values, normalizae the xsec and trials by the number of events
1517  // in the file. This way, we can use it as an approximate event-by-event value
1518  // We do not want to just use the overall value because some of the events may be rejected by various
1519  // event selections, so we only want that ones that were actually use. The easiest way to do so is by
1520  // filling it for each event.
1521  fPythiaTrialsFromFile = trials/nEvents;
1522  // Do __NOT__ divide by nEvents here! The value is already from a TProfile and therefore is already the mean!
1523  fPythiaCrossSectionFromFile = crossSection;
1524 
1525  return true;
1526  }
1527  else {
1528  AliDebugStream(3) << "Unable to open file \"" << pythiaFileName << "\". Will attempt to use values from the hader.";
1529  }
1530 
1531  // Could not open file
1532  return false;
1533 }
1534 
1541 {
1542  if (!fInitializedEmbedding) {
1543  AliError("Chain not initialized before running! Setting up now.");
1544  SetupEmbedding();
1545  }
1546 
1547  // Apply internal event selection
1549  fEmbeddedEventUsed = false;
1550  if (fInternalEventCuts.AcceptEvent(AliAnalysisTaskSE::InputEvent()) == true)
1551  {
1552  fEmbeddedEventUsed = true;
1553  fHistManager.FillTH1("fHistInternalEventCutsStats", "passedEventCuts", 1);
1554 
1555  // The event was accepted by AliEventCuts. Now check for additional cuts.
1556  // Centrality
1557  // NOTE: If the centrality range is the same as AliEventCuts, then simply all will pass
1558  // If a wider centrality range than in AliEventCuts is needed then it must be _entirely_
1559  // configured through manual mode.
1560  if (fCentMin != -999 && fCentMax != -999) {
1561  if (fInternalEventCuts.GetCentrality() < fCentMin || fInternalEventCuts.GetCentrality() > fCentMax) {
1562  fEmbeddedEventUsed = false;
1563  }
1564  else {
1565  fHistManager.FillTH1("fHistInternalEventCutsStats", "centrality", 1);
1566  }
1567  }
1568 
1569  if (fEmbeddedEventUsed) {
1570  // Record all cuts passed
1571  fHistManager.FillTH1("fHistInternalEventCutsStats", "passedAllCuts", 1);
1572  }
1573  }
1574 
1575  // If the internal event was rejected, then record and move on.
1576  if (fEmbeddedEventUsed == false) {
1577  if (fCreateHisto) {
1578  PostData(1, fOutput);
1579  }
1580  return;
1581  }
1582  }
1583 
1584  if (!fInitializedNewFile) {
1585  InitTree();
1586  }
1587 
1588  Bool_t res = GetNextEntry();
1589 
1590  if (!res) {
1591  AliError("Unable to get the event to embed. Nothing will be embedded.");
1592  return;
1593  }
1594 
1595  if (fCreateHisto && fOutput) {
1596  PostData(1, fOutput);
1597  }
1598 }
1599 
1604 {
1605 }
1606 
1612 {
1613  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1614  if (!mgr)
1615  {
1616  AliErrorStream() << "No analysis manager to connect to.\n";
1617  return;
1618  }
1619 
1620  // Remove the dummy task
1621  std::string dummyTaskName = GetName();
1622  dummyTaskName += "_dummyTask";
1623  TObjArray * tasks = mgr->GetTasks();
1624  if (tasks) {
1625  AliAnalysisTaskSE * dummyTask = dynamic_cast<AliAnalysisTaskSE *>(tasks->FindObject(dummyTaskName.c_str()));
1626  if (!dummyTask) {
1627  AliErrorStream() << "Could not remove dummy task \"" << dummyTaskName << "\" from analysis manager! Was it added?\n";
1628  }
1629  // Actually remove the task
1630  tasks->Remove(dummyTask);
1631  AliDebugStream(1) << "Removed dummy task named \"" << dummyTaskName << "\".\n";
1632  }
1633  else {
1634  AliErrorStream() << "Could not retrieve tasks from the analysis manager.\n";
1635  }
1636 }
1637 
1639 {
1641  return &fInternalEventCuts;
1642  }
1643  else {
1644  AliErrorStream() << "Enable manual mode in AliEventCuts (though the embedding helper) to access this object.\n";
1645  }
1646  return nullptr;
1647 }
1648 
1650 {
1651  // Get the pointer to the existing analysis manager via the static access method.
1652  //==============================================================================
1653  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1654  if (!mgr)
1655  {
1656  ::Error("AddTaskEmcalEmbeddingHelper", "No analysis manager to connect to.");
1657  return 0;
1658  }
1659 
1660  // Check the analysis type using the event handlers connected to the analysis manager.
1661  //==============================================================================
1662  AliVEventHandler* handler = mgr->GetInputEventHandler();
1663  if (!handler)
1664  {
1665  ::Error("AddTaskEmcalEmbeddingHelper", "This task requires an input event handler");
1666  return 0;
1667  }
1668 
1669  TString name = "AliAnalysisTaskEmcalEmbeddingHelper";
1670 
1671  AliAnalysisTaskEmcalEmbeddingHelper * mgrTask = static_cast<AliAnalysisTaskEmcalEmbeddingHelper *>(mgr->GetTask(name.Data()));
1672  if (mgrTask) return mgrTask;
1673 
1674  // Create the task that manages
1675  AliAnalysisTaskEmcalEmbeddingHelper * embeddingHelper = new AliAnalysisTaskEmcalEmbeddingHelper(name.Data());
1676 
1677  //-------------------------------------------------------
1678  // Final settings, pass to manager and set the containers
1679  //-------------------------------------------------------
1680 
1681  mgr->AddTask(embeddingHelper);
1682 
1683  // Create containers for input/output
1684  AliAnalysisDataContainer* cInput = mgr->GetCommonInputContainer();
1685 
1686  TString outputContainerName(name);
1687  outputContainerName += "_histos";
1688 
1689  AliAnalysisDataContainer * cOutput = mgr->CreateContainer(outputContainerName.Data(),
1690  TList::Class(),
1691  AliAnalysisManager::kOutputContainer,
1692  Form("%s", AliAnalysisManager::GetCommonFileName()));
1693 
1694  mgr->ConnectInput(embeddingHelper, 0, cInput);
1695  mgr->ConnectOutput(embeddingHelper, 1, cOutput);
1696 
1697  return embeddingHelper;
1698 }
1699 
1701 {
1702  // Get the pointer to the existing analysis manager via the static access method.
1703  //==============================================================================
1704  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1705  if (!mgr)
1706  {
1707  ::Error("ConfigureEmcalEmbeddingHelperOnLEGOTrain", "No analysis manager to connect to.");
1708  return nullptr;
1709  }
1710 
1711  // Retrieve the embedding helper
1712  auto embeddingHelperConst = AliAnalysisTaskEmcalEmbeddingHelper::GetInstance();
1713  // Cast away const-ness on the pointer since the underlying object is not const and we need to be able to modify it.
1714  auto embeddingHelper = const_cast<AliAnalysisTaskEmcalEmbeddingHelper *>(embeddingHelperConst);
1715 
1716  // Fatal if we can't find the task
1717  if (!embeddingHelper) {
1718  AliFatalClass("Could not find embedding helper, Did you remember to create it?");
1719  }
1720 
1721  AliInfoClassStream() << "Found embedding helper to configure.\n";
1722 
1723  // AliAnalysisTaskCfg will require a task to be returned, so we add a dummy task to the analysis manager,
1724  // which will be removed when the user calls Initialize(true) on the embedding helper.
1725  mgr->AddTask(new AliAnalysisTaskSE("AliAnalysisTaskEmcalEmbeddingHelper_dummyTask"));
1726 
1727  return embeddingHelper;
1728 }
1729 
1735 std::string AliAnalysisTaskEmcalEmbeddingHelper::toString(bool includeFileList) const
1736 {
1737  std::stringstream tempSS;
1738 
1739  // General embedding helper information
1740  tempSS << std::boolalpha;
1741  tempSS << GetName() << ": Embedding helper configuration:\n";
1742  tempSS << "Create histos: " << fCreateHisto << "\n";
1743  tempSS << "Pt Hard Bin: " << fPtHardBin << "\n";
1744  tempSS << "N Pt Hard Bins: " << fNPtHardBins << "\n";
1745  tempSS << "File pattern: \"" << fFilePattern << "\"\n";
1746  tempSS << "Input filename: \"" << fInputFilename << "\"\n";
1747  tempSS << "Pythia cross section filename: \"" << fPythiaXSecFilename << "\"\n";
1748  tempSS << "File list filename: \"" << fFileListFilename << "\"\n";
1749  tempSS << "Tree name: " << fTreeName << "\n";
1750  tempSS << "Random event number access: " << fRandomEventNumberAccess << "\n";
1751  tempSS << "Random file access: " << fRandomFileAccess << "\n";
1752  tempSS << "Starting file index: " << fFilenameIndex << "\n";
1753  tempSS << "Number of files to embed: " << fFilenames.size() << "\n";
1754  tempSS << "YAML configuration path: \"" << fConfigurationPath << "\"\n";
1755  tempSS << "Enable internal event selection: " << fUseInternalEventSelection << "\n";
1756  tempSS << "Use manual event cuts mode for internal event selection: " << fUseManualInternalEventCuts << "\n";
1757  if (fCentMin != -999 && fCentMax != -999) {
1758  tempSS << "Internal event selection centrality range: [" << fCentMin << ", " << fCentMax << "]\n";
1759  }
1760  else {
1761  tempSS << "Internal event selection centrality range disabled.\n";
1762  }
1763 
1764  std::bitset<32> triggerMask(fTriggerMask);
1765  tempSS << "\nEmbedded event settings:\n";
1766  tempSS << "Trigger mask (binary): " << triggerMask << "\n";
1767  tempSS << "Reject outliers: " << fMCRejectOutliers << "\n";
1768  tempSS << "Pt hard jet pt rejection factor: " << fPtHardJetPtRejectionFactor << "\n";
1769  tempSS << "Z vertex cut: " << fZVertexCut << "\n";
1770  tempSS << "Max difference between internal and embedded vertex: " << fMaxVertexDist << "\n";
1771 
1772  if (includeFileList) {
1773  tempSS << "\nFiles to embed:\n";
1774  for (auto filename : fFilenames) {
1775  tempSS << "\t" << filename << "\n";
1776  }
1777  }
1778 
1779  return tempSS.str();
1780 }
1781 
1789 std::ostream & AliAnalysisTaskEmcalEmbeddingHelper::Print(std::ostream & in) const {
1790  in << toString();
1791  return in;
1792 }
1793 
1802 std::ostream & operator<<(std::ostream & in, const AliAnalysisTaskEmcalEmbeddingHelper & myTask)
1803 {
1804  std::ostream & result = myTask.Print(in);
1805  return result;
1806 }
1807 
1815 {
1816  std::string temp(opt);
1817  bool includeFileList = false;
1818  if (temp == "FILELIST") {
1819  includeFileList = true;
1820  }
1821  Printf("%s", toString(includeFileList).c_str());
1822 }
Bool_t fRandomEventNumberAccess
If true, it will start embedding from a random entry in the file rather than from the first...
Double_t fZVertexCut
Z vertex cut on embedded event.
const char * filename
Definition: TestFCM.C:1
double Double_t
Definition: External.C:58
Int_t fPtHardBin
ptHard bin for the given pythia production
bool fMCRejectOutliers
If true, MC outliers will be rejected.
double fPythiaCrossSection
! Pythia cross section for the current event (extracted from the pythia header).
Int_t fLowerEntry
! First entry of the current tree to be used for embedding
TString fFileListFilename
Name of the file list containing paths to files to embed.
int fPythiaTrialsFromFile
! Average number of trials extracted from a xsec file.
Bool_t fRandomFileAccess
If true, it will start embedding from a random file in the input files list.
TSystem * gSystem
Int_t fCurrentEntry
! Current entry in the current tree
bool WriteConfiguration(const std::string &filename, const unsigned int i) const
TFile * fExternalFile
! External file used for embedding
Int_t fNPtHardBins
Total number of pt hard bins.
TChain * fChain
! External TChain (tree) containing the events available for embedding
Declaration of class AliAnalysisTaskEmcalEmbeddingHelper.
bool fAutoConfigurePtHardBins
If true, attempt to auto configure pt hard bins. Only works on the LEGO train.
bool fInitializedConfiguration
Notes if the configuration has been initialized.
bool GetProperty(std::vector< std::string > propertyPath, const std::string &propertyName, T &property, const bool requiredProperty) const
bool IsRunInRunlist(const std::string &path) const
AliVHeader * fExternalHeader
! Header of the current external event
AliGenPythiaEventHeader * fPythiaHeader
! Pythia header of the current external event
const AliEventCuts * GetInternalEventCuts() const
Event cuts object for accessing centrality, etc from another task if so inclined. ...
TString fInputFilename
Filename of input root files.
friend std::ostream & operator<<(std::ostream &in, const AliAnalysisTaskEmcalEmbeddingHelper &myTask)
std::vector< std::string > fEmbeddedRunlist
Good runlist for files to embed.
std::string fAutoConfigureBasePath
The base path to the auto configuration (for example, "/alice/cern.ch/user/a/alitrain/") ...
std::string fConfigurationPath
Path to YAML configuration.
static UInt_t DeterminePhysicsSelectionFromYAML(const std::vector< std::string > &selections)
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.
std::string fPythiaXSecFilename
Name of the pythia x sec filename (either "pyxsec.root" or "pyxsec_hists.root")
unsigned int UInt_t
Definition: External.C:33
THashList * GetListOfHistograms() const
Get the list of histograms.
Definition: THistManager.h:671
float Float_t
Definition: External.C:68
bool fInitializedEmbedding
! Notes where the TChain has been initialized for embedding
bool fUseInternalEventSelection
If true, apply internal event selection though AliEventCuts.
std::string toString(bool includeFileList=false) const
Implementation of task to embed external events.
UInt_t fMaxNumberOfFiles
! Max number of files that are in the TChain
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
Create a new TH1 within the container.
void FillProfile(const char *name, double x, double y, double weight=1.)
Double_t nEvents
plot quality messages
Double_t fMaxVertexDist
Max distance between Z vertex of internal and embedded event.
AliVEvent * fExternalEvent
! Current external event available for embedding
std::vector< std::string > fPythiaCrossSectionFilenames
Paths to the pythia xsection files.
void FillTH1(const char *hname, double x, double weight=1., Option_t *opt="")
Fill a 1D histogram within the container.
TString fTreeName
Name of the ESD/AOD tree where the events are to be found.
Int_t fOffset
! Offset from fLowerEntry where the loop over the tree should start
bool fCreateHisto
If true, create QA histograms.
bool fEmbeddedEventUsed
! If true, the internal event was selected, so the embedded event is used. Defaults to true so other ...
Enhanced TList-derived class that implements correct merging for pt_hard binned production.
Definition: AliEmcalList.h:25
bool fWrappedAroundTree
! Notes whether we have wrapped around the tree, which is important if the offset into the tree is no...
Double_t fPtHardJetPtRejectionFactor
Factor which the pt hard bin is multiplied by to compare against pythia header jets pt...
PWG::Tools::AliYAMLConfiguration fYAMLConfig
Hanldes configuration from YAML.
int AddEmptyConfiguration(const std::string &configurationName)
Add YAML configuration at configurationFilename to available configurations.
double fPythiaPtHard
! Pt hard of the current event (extracted from the pythia header).
static AliAnalysisTaskEmcalEmbeddingHelper * AddTaskEmcalEmbeddingHelper()
THistManager fHistManager
Manages access to all histograms.
bool fInitializedNewFile
! Notes where the entry indices have been initialized for a new tree in the chain ...
static AliAnalysisTaskEmcalEmbeddingHelper * fgInstance
! Global instance of this class
UInt_t fFileNumber
! File number corresponding to the current tree
bool fUseManualInternalEventCuts
If true, manual event cuts mode will be used for AliEventCuts.
bool IsFileAccessible(std::string filename)
YAML configuration class for AliPhysics.
bool WriteProperty(std::string propertyName, T &property, std::string configurationName="")
std::string RemoveTrailingSlashes(std::string filename) const
int AddConfiguration(std::string configurationFilename, std::string configurationName="")
double fCentMax
Maximum centrality for internal event selection.
const char Option_t
Definition: External.C:48
AliEmcalList * fOutput
! List which owns the output histograms to be saved
AliEventCuts fInternalEventCuts
If enabled, Handles internal event selection.
bool Bool_t
Definition: External.C:53
std::string fAutoConfigureIdentifier
How the auto configuration YAML file should be identified. (for example, "rehlersTrain") ...
Int_t fUpperEntry
! Last entry of the current tree to be used for embedding
Int_t fFilenameIndex
Index of vector containing paths to files to embed.
double fPythiaCrossSectionFromFile
! Average pythia cross section extracted from a xsec file.
std::vector< std::string > fFilenames
Paths to the files to embed.
static AliAnalysisTaskEmcalEmbeddingHelper * ConfigureEmcalEmbeddingHelperOnLEGOTrain()
std::string fAutoConfigureTrainTypePath
The path associated with the train type (for example, "PWGJE/Jets_EMC_PbPb/")
TString fFilePattern
File pattern to select AliEn files using alien_find.
int fPythiaTrials
! Number of pythia trials for the current event (extracted from the pythia header).
double fCentMin
Minimum centrality for internal event selection.
TList * OpenFile(const char *fname)
Definition: DrawAnaELoss.C:65
static const AliAnalysisTaskEmcalEmbeddingHelper * GetInstance()
std::string ConstructFullPythiaXSecFilename(std::string inputFilename, const std::string &pythiaFilename, bool testIfExists) const