AliPhysics  5eaf189 (5eaf189)
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  fPrintTimingInfoToLog(false),
168  fTimer()
169 {
170  if (fgInstance != nullptr) {
171  AliError("An instance of AliAnalysisTaskEmcalEmbeddingHelper already exists: it will be deleted!!!");
172  delete fgInstance;
173  }
174 
175  fgInstance = this;
176 }
177 
184  AliAnalysisTaskSE(name),
185  fTriggerMask(0),
186  fMCRejectOutliers(false),
188  fZVertexCut(10),
189  fMaxVertexDist(999),
191  fInitializedNewFile(false),
192  fInitializedEmbedding(false),
193  fWrappedAroundTree(false),
194  fTreeName("aodTree"),
195  fNPtHardBins(1),
196  fPtHardBin(-1),
197  fRandomEventNumberAccess(kFALSE),
198  fRandomFileAccess(kTRUE),
199  fCreateHisto(true),
200  fYAMLConfig(),
204  fEmbeddedEventUsed(true),
205  fCentMin(-999),
206  fCentMax(-999),
208  fAutoConfigureBasePath("alien:///alice/cern.ch/user/a/alitrain/"),
209  fAutoConfigureTrainTypePath("PWGJE/Jets_EMC_PbPb/"),
210  fAutoConfigureIdentifier("autoConfigIdentifier"),
211  fFilePattern(""),
212  fInputFilename(""),
213  fFileListFilename(""),
214  fFilenameIndex(-1),
215  fFilenames(),
216  fConfigurationPath(""),
220  fChain(nullptr),
221  fCurrentEntry(0),
222  fLowerEntry(0),
223  fUpperEntry(0),
224  fOffset(0),
226  fFileNumber(0),
227  fHistManager(name),
228  fOutput(nullptr),
232  fPythiaTrials(0),
236  fPythiaPtHard(0.),
237  fPrintTimingInfoToLog(false),
238  fTimer()
239 {
240  if (fgInstance != 0) {
241  AliError("An instance of AliAnalysisTaskEmcalEmbeddingHelper already exists: it will be deleted!!!");
242  delete fgInstance;
243  }
244 
245  fgInstance = this;
246 
247  if (fCreateHisto) {
248  DefineOutput(1, AliEmcalList::Class());
249  }
250 }
251 
258 {
259  if (fgInstance == this) fgInstance = nullptr;
260  if (fExternalEvent) delete fExternalEvent;
261  if (fExternalFile) {
262  fExternalFile->Close();
263  delete fExternalFile;
264  }
265 }
266 
268 {
269  // Initialize %YAML configuration, if one is given
270  bool initializedYAML = InitializeYamlConfig();
271 
273 
274  // Get file list
275  bool result = GetFilenames();
276 
277  if (result && initializedYAML) {
279  }
280 
281  if (removeDummyTask == true) {
282  RemoveDummyTask();
283  }
284 
285  // Initialize the YAML config object for streaming
287 
288  // Print the results of the initialization
289  // Print outside of the ALICE Log system to ensure that it is always available!
290  std::cout << *this;
291 
292  return result;
293 }
294 
302 {
303  // Following the variable blocks defined in the header
304  // Embedded event properties
305  std::vector<std::string> physicsSelection;
306  bool res = fYAMLConfig.GetProperty("embeddedEventPhysicsSelection", physicsSelection, false);
307  if (res) {
309  }
310  res = fYAMLConfig.GetProperty("enableMCOutlierRejection", fMCRejectOutliers, false);
311  res = fYAMLConfig.GetProperty("ptHardJetPtRejectionFactor", fPtHardJetPtRejectionFactor, false);
312  res = fYAMLConfig.GetProperty("embeddedEventZVertexCut", fZVertexCut, false);
313  res = fYAMLConfig.GetProperty("maxVertexDifferenceDistance", fMaxVertexDist, false);
314 
315  // Embedding helper properties
316  res = fYAMLConfig.GetProperty("treeName", fTreeName, false);
317  res = fYAMLConfig.GetProperty("nPtHardBins", fNPtHardBins, false);
318  res = fYAMLConfig.GetProperty("ptHardBin", fPtHardBin, false);
319  res = fYAMLConfig.GetProperty("randomEventNumberAccess", fRandomEventNumberAccess, false);
320  res = fYAMLConfig.GetProperty("randomFileAccess", fRandomFileAccess, false);
321  res = fYAMLConfig.GetProperty("createHisto", fCreateHisto, false);
322  res = fYAMLConfig.GetProperty("printTimingInfoInLog", fPrintTimingInfoToLog, false);
323  // More general embedding helper properties
324  res = fYAMLConfig.GetProperty("filePattern", fFilePattern, false);
325  res = fYAMLConfig.GetProperty("inputFilename", fInputFilename, false);
326  res = fYAMLConfig.GetProperty("fileListFilename", fFileListFilename, false);
327  res = fYAMLConfig.GetProperty("filenameIndex", fFilenameIndex, false);
328  // Configuration path makes no sense, as we are already using the %YAML configuration
329  res = fYAMLConfig.GetProperty("runlist", fEmbeddedRunlist, false);
330  // Generally should not be set
331  res = fYAMLConfig.GetProperty("filenames", fFilenames, false);
332  res = fYAMLConfig.GetProperty("fPythiaCrossSectionFilenames", fPythiaCrossSectionFilenames, false);
333 
334  // Internal event selection properties
335  // NOTE: Need to define the base name here so that the property path is not ambiguous (due to otherwise only being `const char *`)
336  std::string baseName = "internalEventSelection";
337  res = fYAMLConfig.GetProperty({baseName, "enabled"}, fUseInternalEventSelection, false);
338  res = fYAMLConfig.GetProperty({baseName, "useManualCuts"}, fUseManualInternalEventCuts, false);
339  // Centrality
340  std::vector <double> centralityRange;
341  res = fYAMLConfig.GetProperty({baseName, "centralityRange"}, centralityRange, false);
342  if (res) {
343  if (centralityRange.size() != 2) {
344  AliErrorStream() << "Passed centrality range with " << centralityRange.size() << " entries, but 2 values are required. Ignoring values.\n";
345  }
346  else {
347  AliDebugStream(1) << "Setting internal event centrality range to [" << centralityRange.at(0) << ", " << centralityRange.at(1) << "]\n";
348  fCentMin = centralityRange.at(0);
349  fCentMax = centralityRange.at(1);
350  }
351  }
352  // Physics selection
353  res = fYAMLConfig.GetProperty({baseName, "physicsSelection"}, physicsSelection, false);
354  if (res) {
355  fOfflineTriggerMask = AliEmcalContainerUtils::DeterminePhysicsSelectionFromYAML(physicsSelection);
356  }
357 
358  // Auto configure pt hard properties
359  res = fYAMLConfig.GetProperty("autoConfigurePtHardBins", fAutoConfigurePtHardBins, false);
360  res = fYAMLConfig.GetProperty("autoConfigureBasePath", fAutoConfigureBasePath, false);
361  res = fYAMLConfig.GetProperty("autoConfigureTrainTypePath", fAutoConfigureTrainTypePath, false);
362  res = fYAMLConfig.GetProperty("autoConfigureIdentifier", fAutoConfigureIdentifier, false);
363 }
364 
389 {
390  // Determine the pattern filename if not yet set
391  if (fInputFilename == "") {
392  if (fTreeName == "aodTree") {
393  fInputFilename = "AliAOD.root";
394  }
395  else if (fTreeName == "esdTree") {
396  fInputFilename = "AliESDs.root";
397  }
398  else {
399  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()));
400  }
401  }
402 
403  // Retrieve filenames if we don't have them yet.
404  if (fFilenames.size() == 0)
405  {
406  // Handle pt hard bin auto configuration
408  {
409  if (fPtHardBin > 0) {
410  AliFatal("Requested both pt hard bin auto configuration and selected a non-zero pt hard bin. These are incompatible options. Please check your configuration.");
411  }
412  bool success = AutoConfigurePtHardBins();
413  if (success == false) {
414  AliFatal("Pt hard bin auto configuration requested, but it failed. Please check the logs.\n");
415  }
416  }
417 
418  // Handle if fPtHardBin is set
419  // This will require formatting the file pattern in the proper way to support these substitutions
420  if (fPtHardBin != -1 && fFilePattern != "") {
421  fFilePattern = TString::Format(fFilePattern, fPtHardBin);
422  }
423 
424  // Setup AliEn access if needed
425  if (fFilePattern.Contains("alien://") || fFileListFilename.Contains("alien://")) {
427  }
428 
429  // Retrieve AliEn filenames directly from AliEn
430  bool usedFilePattern = false;
431  if (fFilePattern.Contains("alien://")) {
432  usedFilePattern = true;
433  AliDebug(2, TString::Format("Trying to retrieve file list from AliEn with pattern file %s...", fFilePattern.Data()));
434 
435  // Create a temporary filename based on a UUID to make sure that it doesn't overwrite any files
436  if (fFileListFilename == "") {
438  }
439 
440  // The query command cannot handle "alien://" in the file pattern, so we need to remove it for the command
441  TString filePattern = fFilePattern;
442  filePattern.ReplaceAll("alien://", "");
443 
444  // Execute the grid query to get the filenames
445  AliDebug(2, TString::Format("Trying to retrieve file list from AliEn with pattern \"%s\" and input filename \"%s\"", filePattern.Data(), fInputFilename.Data()));
446  auto result = gGrid->Query(filePattern.Data(), fInputFilename.Data());
447 
448  if (result) {
449 
450  // Loop over the result to store it in the fileList file
451  std::ofstream outFile(fFileListFilename);
452  for (int i = 0; i < result->GetEntries(); i++)
453  {
454  TString path = result->GetKey(i, "turl");
455  // "turl" corresponds to the full AliEn url
456 
457  // If a runlist is specified for good embedded runs, only include the file if it is in this runlist
458  if (IsRunInRunlist(path.Data())) {
459  outFile << path << "\n";
460  }
461  }
462  outFile.close();
463  }
464  else {
465  AliErrorStream() << "Failed to run grid query\n";
466  return false;
467  }
468  }
469 
470  // Handle a filelist on AliEn
471  if (fFileListFilename.Contains("alien://")) {
472  // Check if we already used the file pattern
473  if (usedFilePattern) {
474  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";
475  }
476 
477  // Determine the local filename and copy file to local directory
478  std::string alienFilename = fFileListFilename.Data();
479  fFileListFilename = gSystem->BaseName(alienFilename.c_str());
481 
482  TFile::Cp(alienFilename.c_str(), fFileListFilename.Data());
483  }
484 
485  std::ifstream inputFile(fFileListFilename);
486 
487  // Copy available files into the filenames vector
488  // From:: https://stackoverflow.com/a/8365247
489  std::copy(std::istream_iterator<std::string>(inputFile),
490  std::istream_iterator<std::string>(),
491  std::back_inserter(fFilenames));
492 
493  inputFile.close();
494  }
495 
496  if (fFilenames.size() == 0) {
497  AliFatal(TString::Format("Filenames from pattern \"%s\" and file list \"%s\" yielded an empty list!", fFilePattern.Data(), fFileListFilename.Data()));
498  }
499 
500  // Add "#" to files in there are any zip files
501  // It is require to open the proper root file within the zip
502  for (auto filename : fFilenames)
503  {
504  if (filename.find(".zip") != std::string::npos && filename.find("#") == std::string::npos) {
505  filename += "#";
507  if (fTreeName == "aodTree") {
508  filename += "#AliAOD.root";
509  }
510  else if (fTreeName == "esdTree") {
511  filename += "#AliESDs.root";
512  }
513  else {
514  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()));
515  return false;
516  }
517  }
518  }
519 
520  // Determine whether AliEn is needed
521  // It is possible that this has not been determined up to this point
522  for (auto filename : fFilenames)
523  {
524  if (filename.find("alien://") != std::string::npos) {
526  // No point in continuing to search once we know that it is needed
527  break;
528  }
529  }
530 
531  // Check if each filenames exists. If they do not exist, then remove them for fFilenames
532  unsigned int initialSize = fFilenames.size();
533  // NOTE: We invert the result of IsFileAccessible because we should return true for files that should be _removed_ (ie are inaccessible)
534  fFilenames.erase(std::remove_if(fFilenames.begin(), fFilenames.end(), [](const std::string & filename) {return (::IsFileAccessible(filename) == false);} ), fFilenames.end());
535 
536  AliInfoStream() << "Found " << fFilenames.size() << " files to embed (" << (initialSize - fFilenames.size()) << " filename(s) inaccessible or invalid)\n";
537 
538  // Determine pythia filename
540 
541  return true;
542 }
543 
549 {
550  // Get the initial filename. Use the first entry as a proxy for other input files
551  std::string externalEventFilename = "";
552  if (fFilenames.size() > 0) {
553  externalEventFilename = fFilenames.at(0);
554  }
555  else {
556  return;
557  }
558 
559  std::vector <std::string> pythiaBaseFilenames = {"pyxsec.root", "pyxsec_hists.root"};
560  AliInfoStream() << "Attempting to determine pythia cross section filename. It can be normal to see some TFile::Init() errors!\n";
561  std::string pythiaXSecFilename = "";
562  for (auto & name : pythiaBaseFilenames) {
563  pythiaXSecFilename = ConstructFullPythiaXSecFilename(externalEventFilename, name, true);
564  if (pythiaXSecFilename != "") {
565  AliDebugStream(4) << "Found pythia cross section filename \"" << name << "\"\n";
566  fPythiaXSecFilename = name;
567  break;
568  }
569  }
570 
571  if (fPythiaXSecFilename == "") {
572  // Failed entirely - just give up on this
573  // We will use an empty filename as a proxy for whether the file has been found (empty is equivalent to not found)
574  AliErrorStream() << "Failed to find pythia x sec file! Continuing with only the pythia header!\n";
575  }
576  else {
577  AliInfoStream() << "Found pythia cross section file \"" << fPythiaXSecFilename << "\".\n";
578  }
579 }
580 
581 
589 bool AliAnalysisTaskEmcalEmbeddingHelper::IsRunInRunlist(const std::string & path) const
590 {
591  if (fEmbeddedRunlist.size() == 0) {
592  return true;
593  }
594 
595  for (auto run : fEmbeddedRunlist) {
596  if (path.find(run) != std::string::npos) {
597  return true;
598  }
599  }
600  return false;
601 }
602 
609 {
610  if (fConfigurationPath == "") {
611  AliInfo("No Embedding YAML configuration was provided");
612  }
613  else {
614  AliInfoStream() << "Embedding YAML configuration was provided: \"" << fConfigurationPath << "\".\n";
615 
616  int addedConfig = fYAMLConfig.AddConfiguration(fConfigurationPath, "yamlConfig");
617  if (addedConfig < 0) {
618  AliError("YAML Configuration not found!");
619  return false;
620  }
621  }
622 
623  return true;
624 }
625 
642 {
643  bool returnValue = false;
644 
645  AliInfoStream() << "Attempting to auto configure pt hard bins.\n";
646  // %YAML configuration containing pt hard bin to train number mapping
648 
649  // Handle AliEn explicitly here since the default base path contains "alien://"
650  if (fAutoConfigureBasePath.find("alien://") != std::string::npos && !gGrid) {
652  }
653 
654  // Get train ID
655  // Need to get the char * directly because it may be null.
656  const char * trainNumberStr = gSystem->Getenv("TRAIN_RUN_ID");
657  std::stringstream trainNumberSS;
658  if (trainNumberStr) {
659  trainNumberSS << trainNumberStr;
660  }
661  if (trainNumberSS.str() == "") {
662  AliFatal("Cannot retrieve train ID.");
663  }
664  // Extract train number from the string
665  int trainNumber;
666  trainNumberSS >> trainNumber;
667 
668  // Determine the file path
670  filename += "/";
672  filename += "/";
674  // Add ".yaml" if it is not already there
675  std::string yamlExtension = ".yaml";
676  if (filename.find(yamlExtension) == std::string::npos) {
677  filename += yamlExtension;
678  }
679 
680  // Check if file exists
681  if (gSystem->AccessPathName(filename.c_str())) {
682  // File _does not_ exist
683  AliInfoStream() << "Train pt hard bin configuration file not available, so creating a new empty configuration named \"" << fAutoConfigureIdentifier << "\".\n";
684  // Use an empty configuration
686  }
687  else {
688  AliInfoStream() << "Opening configuration located at \"" << filename << "\".\n";
689  // Use the existing configuration
691  }
692 
693  // Look for each pt hard bin, and then retrieve the corresponding train number
694  // Once an open pt hard bin is found, add the current train number
695  int tempTrainNumber = -1;
696  bool getPropertyReturnValue = false;
697  std::stringstream propertyName;
698  for (int ptHardBin = 1; ptHardBin <= fNPtHardBins; ptHardBin++)
699  {
700  propertyName.str("");
701  propertyName << ptHardBin;
702  getPropertyReturnValue = config.GetProperty(propertyName.str(), tempTrainNumber, false);
703  if (getPropertyReturnValue != true) {
704  AliInfoStream() << "Train " << trainNumber << " will use pt hard bin " << ptHardBin << ".\n";
705  // We have determine our pt hard bin!
706  fPtHardBin = ptHardBin;
707 
708  // Write the train number back out to the %YAML configuration and save it
709  config.WriteProperty(propertyName.str(), trainNumber, fAutoConfigureIdentifier);
711 
712  // NOTE: Cannot clean up the yaml file on the last pt hard bin because the train can be launched
713  // multiple times due to tests, etc. Therefore, we have to accept that we are leaving around used
714  // yaml config files.
715 
716  // We are done - continue on.
717  returnValue = true;
718  break;
719  }
720  else {
721  AliDebugStream(2) << "Found pt hard bin " << ptHardBin << " corresponding to train number " << trainNumber << ".\n";
722  // If train was already allocated (say, by a test train), then use that pt hard bin
723  if (tempTrainNumber == trainNumber) {
724  AliInfoStream() << "Train run number " << trainNumber << " was already found assigned to pt hard bin " << ptHardBin << ". That pt hard bin will be used.\n";
725  fPtHardBin = ptHardBin;
726 
727  // We are done - continue on.
728  returnValue = true;
729  break;
730  }
731  // Otherwise, nothing to be done.
732  }
733  }
734 
735  return returnValue;
736 }
737 
746 {
747  std::string tempStr = "";
748  if (fFileListFilename == "") {
749  tempStr = "fileList";
750  }
751  TUUID tempUUID;
752  tempStr += ".";
753  tempStr += tempUUID.AsString();
754  tempStr += ".txt";
755 
756  return tempStr;
757 }
758 
767 {
768  while (filename.rbegin() != filename.rend() && *(filename.rbegin()) == '/') {
769  filename.pop_back();
770  }
771 
772  return filename;
773 }
774 
780 {
781  // This determines which file is added first to the TChain, thus determining the order of processing
782  // Random file access. Only do this if the user has no set the filename index and request random file access
783  if (fFilenameIndex == -1 && fRandomFileAccess) {
784  // Floor ensures that we it doesn't overflow
785  TRandom3 rand(0);
786  fFilenameIndex = TMath::FloorNint(rand.Rndm()*fFilenames.size());
787  // +1 to account for the fact that the filenames vector is 0 indexed.
788  AliInfo(TString::Format("Starting with random file number %i!", fFilenameIndex+1));
789  }
790  // If not random file access, then start from the beginning
791  if (fFilenameIndex < 0 || static_cast<UInt_t>(fFilenameIndex) >= fFilenames.size()) {
792  // Skip notifying on -1 since it will likely be set there due to constructor.
793  if (fFilenameIndex != -1) {
794  AliWarning(TString::Format("File index %i out of range from 0 to %lu! Resetting to 0!", fFilenameIndex, fFilenames.size()));
795  }
796  fFilenameIndex = 0;
797  }
798 
799  // +1 to account for the fact that the filenames vector is 0 indexed.
800  AliInfo(TString::Format("Starting with file number %i out of %lu", fFilenameIndex+1, fFilenames.size()));
801 }
802 
812 {
813  Int_t attempts = -1;
814 
815  do {
816  // Reset to start of tree
817  if (fCurrentEntry == fUpperEntry) {
819  fWrappedAroundTree = true;
820  }
821 
823  // Continue with GetEntry as normal
824  }
825  else {
826  // NOTE: On transition from one file to the next, this calls the next entry that would be expected.
827  // However, if it is for the last file, it tries to GetEntry() of one entry past the end of the last file.
828  // Normally, this would be a problem, however GetEntry() just doesn't fill the fields of an invalid index
829  // instead of throwing an error. So "invalid values" are filled for a file that doesn't exist, but then
830  // they are immediately replaced by the lines below that reset the access values and re-init the tree.
831  // The benefit of this approach is it simplies file counting (we don't need to carefully increment here
832  // and in InitTree()) and preserves the desired behavior when we are not at the last file.
833  InitTree();
834  }
835 
836  // Load current event
837  // Can be a simple less than, because fFileNumber counts from 0.
839  fChain->GetEntry(fCurrentEntry);
840  }
841  else {
842  AliError("====================================================================================================");
843  AliError("== No more files available to embed from the TChain! Restarting from the beginning of the TChain! ==");
844  AliError("== Be careful to check that this is the desired action! ==");
845  AliError("====================================================================================================");
846 
847  // Reset the relevant access values
848  // fCurrentEntry and fLowerEntry are automatically reset in InitTree()
849  fFileNumber = 0;
850  fUpperEntry = 0;
851 
852  // Re-init back to the start
853  InitTree();
854 
855  // Access the relevant entry
856  // We are certain that fFileNumber is less than fMaxNumberOfFiles, so we are resetting to start
857  fChain->GetEntry(fCurrentEntry);
858  }
859  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));
860 
861  // Set relevant event properties
863 
864  // Increment current entry
865  fCurrentEntry++;
866 
867  // Provide a check for number of attempts
868  attempts++;
869  if (attempts == 1000)
870  AliWarning("After 1000 attempts no event has been accepted by the event selection (trigger, centrality...)!");
871 
872  // Record event properties
873  if (fCreateHisto) {
875  }
876 
877  } while (!IsEventSelected());
878 
879  if (fCreateHisto) {
880  fHistManager.FillTH1("fHistEventCount", "Accepted");
881  fHistManager.FillTH1("fHistEmbeddedEventsAttempted", attempts);
882  }
883 
884  if (!fChain) return kFALSE;
885 
886  return kTRUE;
887 }
888 
894 {
895  AliDebug(4, "Set event properties");
896  fExternalHeader = fExternalEvent->GetHeader();
897 
898  // Handle pythia header if AOD
899  AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(fExternalEvent->FindListObject(AliAODMCHeader::StdBranchName()));
900  if (aodMCH) {
901  for (UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++) {
902  fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
903  if (fPythiaHeader) break;
904  }
905  }
906 
907  if (fPythiaHeader)
908  {
909  fPythiaCrossSection = fPythiaHeader->GetXsection();
910  fPythiaTrials = fPythiaHeader->Trials();
911  fPythiaPtHard = fPythiaHeader->GetPtHard();
912  // It is identically zero if the cross section is not available
913  if (fPythiaCrossSection == 0.) {
914  AliDebugStream(4) << "Taking the pythia cross section avg from the xsec file.\n";
916  }
917  // It is identically zero if the number of trials is not available
918  if (fPythiaTrials == 0.) {
919  AliDebugStream(4) << "Taking the pythia trials avg from the xsec file.\n";
921  }
922  // Pt hard is inherently event-by-event and cannot by taken as a avg quantity.
923 
924  AliDebugStream(4) << "Pythia header is defined!\n";
925  AliDebugStream(4) << "fPythiaCrossSection: " << fPythiaCrossSection << "\n";
926  }
927 }
928 
933 {
934  // Fill trials, xsec, pt hard
935  fHistManager.FillTH1("fHistTrials", fPtHardBin, fPythiaTrials);
937  fHistManager.FillTH1("fHistPtHard", fPythiaPtHard);
938 }
939 
946 {
948  return kTRUE;
949  }
950 
951  if (fCreateHisto) {
952  // Keep count of number of rejected events
953  fHistManager.FillTH1("fHistEventCount", "Rejected");
954  }
955 
956  return kFALSE;
957 }
958 
965 {
966  // Check if pt hard bin is 0, indicating a problem with the event or the grid.
967  // In such a case, the event should be rejected.
968  // This condition should only be applied if we have a valid pythia header.
969  // (pt hard should still be set even if the production wasn't done in pt hard bins).
970  if (fPythiaPtHard == 0. && fPythiaHeader) {
971  AliDebugStream(3) << "Event rejected due to pt hard = 0, indicating a problem with the external event.\n";
972  if (fCreateHisto) {
973  fHistManager.FillTH1("fHistEmbeddedEventRejection", "PtHardIs0", 1);
974  }
975  return kFALSE;
976  }
977 
978  // Physics selection
979  if (fTriggerMask != 0) {
980  UInt_t res = 0;
981  const AliESDEvent *eev = dynamic_cast<const AliESDEvent*>(fExternalEvent);
982  if (eev) {
983  AliFatal("Event selection is not implemented for embedding ESDs.");
984  // Unfortunately, the normal method of retrieving the trigger mask (commented out below) doesn't work for the embedded event since we don't
985  // 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
986  // probably best to avoid it if possible.
987  //
988  // Suggestions are welcome here!
989  //res = (dynamic_cast<AliInputEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
990  } else {
991  const AliAODEvent *aev = dynamic_cast<const AliAODEvent*>(fExternalEvent);
992  if (aev) {
993  res = (dynamic_cast<AliVAODHeader*>(aev->GetHeader()))->GetOfflineTrigger();
994  }
995  }
996 
997  if ((res & fTriggerMask) == 0) {
998  AliDebug(3, Form("Event rejected due to physics selection. Event trigger mask: %d, trigger mask selection: %d.",
999  res, fTriggerMask));
1000  if (fCreateHisto) {
1001  fHistManager.FillTH1("fHistEmbeddedEventRejection", "PhysSel", 1);
1002  }
1003  return kFALSE;
1004  }
1005  }
1006 
1007  // Vertex selection
1008  Double_t externalVertex[3]={0};
1009  Double_t inputVertex[3]={0};
1010  const AliVVertex *externalVert = fExternalEvent->GetPrimaryVertex();
1011  const AliVVertex *inputVert = AliAnalysisTaskSE::InputEvent()->GetPrimaryVertex();
1012  if (externalVert && inputVert) {
1013  externalVert->GetXYZ(externalVertex);
1014  inputVert->GetXYZ(inputVertex);
1015 
1016  if (TMath::Abs(externalVertex[2]) > fZVertexCut) {
1017  AliDebug(3, Form("Event rejected due to Z vertex selection. Event Z vertex: %f, Z vertex cut: %f",
1018  externalVertex[2], fZVertexCut));
1019  if (fCreateHisto) {
1020  fHistManager.FillTH1("fHistEmbeddedEventRejection", "Vz", 1);
1021  }
1022  return kFALSE;
1023  }
1024  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]));
1025  if (dist > fMaxVertexDist) {
1026  AliDebug(3, Form("Event rejected because the distance between the current and embedded vertices is > %f. "
1027  "Current event vertex (%f, %f, %f), embedded event vertex (%f, %f, %f). Distance = %f",
1028  fMaxVertexDist, inputVertex[0], inputVertex[1], inputVertex[2], externalVertex[0], externalVertex[1], externalVertex[2], dist));
1029  if (fCreateHisto) {
1030  fHistManager.FillTH1("fHistEmbeddedEventRejection", "VertexDist", 1);
1031  }
1032  return kFALSE;
1033  }
1034  }
1035 
1036  // Check for pt hard bin outliers
1038  {
1039  // Pythia jet / pT-hard > factor
1040  // This corresponds to "condition 1" in AliAnalysisTaskEmcal
1041  // NOTE: The other "conditions" defined there are not really suitable to define here, since they
1042  // depend on the input objects of the event
1043  if (fPtHardJetPtRejectionFactor > 0.) {
1044  TLorentzVector jet;
1045 
1046  Int_t nTriggerJets = fPythiaHeader->NTriggerJets();
1047 
1048  AliDebugStream(4) << "Pythia Njets: " << nTriggerJets << ", pT Hard: " << fPythiaPtHard << "\n";
1049 
1050  Float_t tmpjet[]={0,0,0,0};
1051  for (Int_t iJet = 0; iJet< nTriggerJets; iJet++) {
1052  fPythiaHeader->TriggerJet(iJet, tmpjet);
1053 
1054  jet.SetPxPyPzE(tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3]);
1055 
1056  AliDebugStream(5) << "Pythia jet " << iJet << ", pycell jet pT: " << jet.Pt() << "\n";
1057 
1058  //Compare jet pT and pt Hard
1059  if (jet.Pt() > fPtHardJetPtRejectionFactor * fPythiaPtHard) {
1060  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";
1061  fHistManager.FillTH1("fHistEmbeddedEventRejection", "MCOutlier", 1);
1062  return kFALSE;
1063  }
1064  }
1065  }
1066  }
1067 
1068  return kTRUE;
1069 }
1070 
1077 {
1078  if (!fChain) return kFALSE;
1079 
1080  if (!fExternalEvent) {
1081  if (fTreeName == "aodTree") {
1082  fExternalEvent = new AliAODEvent();
1083  }
1084  else if (fTreeName == "esdTree") {
1085  fExternalEvent = new AliESDEvent();
1086  }
1087  else {
1088  AliError(Form("Tree name %s not recognized!", fTreeName.Data()));
1089  return kFALSE;
1090  }
1091  }
1092 
1093  fExternalEvent->ReadFromTree(fChain, fTreeName);
1094 
1095  return kTRUE;
1096 }
1097 
1104 {
1105  SetupEmbedding();
1106 
1107  // Reinitialize the YAML config after it was streamed so that it can be used properly.
1109 
1110  // Setup AliEventCuts
1112  {
1113  AliDebugStream(1) << "Configuring AliEventCuts for internal event selection.\n";
1114  // Handle manual cuts
1116  fInternalEventCuts.SetManualMode();
1117  // Implement these cuts by retrieving the event cuts object and setting them manually.
1118  }
1119 
1120  // Trigger selection
1121  bool useEventCutsAutomaticTriggerSelection = false;
1122  bool res = fYAMLConfig.GetProperty(std::vector<std::string>({"internalEventSelection", "useEventCutsAutomaticTriggerSelection"}), useEventCutsAutomaticTriggerSelection, false);
1123  if (res && useEventCutsAutomaticTriggerSelection) {
1124  // Use the autmoatic selection. Nothing to be done.
1125  AliDebugStream(1) << "Using the automatic trigger selection from AliEventCuts.\n";
1126  }
1127  else {
1128  // Use the cuts selected by SelectCollisionCandidates()
1129  AliDebugStream(1) << "Using the trigger selection specified with SelectCollisionCandidates()\n.";
1130  fInternalEventCuts.OverrideAutomaticTriggerSelection(fOfflineTriggerMask);
1131  }
1132  }
1133 
1134  // Set up timer for logging purposes
1135  if (fPrintTimingInfoToLog) {
1136  fTimer = TStopwatch();
1137  }
1138 
1139  if (!fCreateHisto) {
1140  return;
1141  }
1142 
1143  // Create output list
1144  OpenFile(1);
1145  fOutput = new AliEmcalList();
1146  fOutput->SetOwner();
1147 
1148  // Get the histograms from AliEventCuts
1150  // This list will be owned by fOutput, so it won't be leaked.
1151  TList * eventCutsOutput = new TList();
1152  eventCutsOutput->SetOwner(kTRUE);
1153  eventCutsOutput->SetName("EventCuts");
1154 
1155  // Add the event cuts to the output
1156  fInternalEventCuts.AddQAplotsToList(eventCutsOutput);
1157  fOutput->Add(eventCutsOutput);
1158  }
1159 
1160  // Create histograms
1161  TString histName;
1162  TString histTitle;
1163 
1164  // Cross section
1165  histName = "fHistXsection";
1166  histTitle = "Pythia Cross Section;p_{T} hard bin; XSection";
1167  fHistManager.CreateTProfile(histName, histTitle, fNPtHardBins, 0, fNPtHardBins);
1168 
1169  // Trials
1170  histName = "fHistTrials";
1171  histTitle = "Number of Pythia Trials;p_{T} hard bin;Trials";
1172  fHistManager.CreateTH1(histName, histTitle, fNPtHardBins, 0, fNPtHardBins);
1173 
1174  // Pt hard spectra
1175  histName = "fHistPtHard";
1176  histTitle = "p_{T} Hard Spectra;p_{T} hard;Counts";
1177  fHistManager.CreateTH1(histName, histTitle, 500, 0, 1000);
1178 
1179  // Count of accepted and rejected events
1180  histName = "fHistEventCount";
1181  histTitle = "Event count;Result;Count";
1182  auto histEventCount = fHistManager.CreateTH1(histName, histTitle, 2, 0, 2);
1183  histEventCount->GetXaxis()->SetBinLabel(1,"Accepted");
1184  histEventCount->GetXaxis()->SetBinLabel(2,"Rejected");
1185 
1186  // Event rejection reason
1187  histName = "fHistEmbeddedEventRejection";
1188  histTitle = "Reasons to reject embedded event";
1189  std::vector<std::string> binLabels = {"PhysSel", "MCOutlier", "Vz", "VertexDist", "PtHardIs0"};
1190  auto histEmbeddedEventRejection = fHistManager.CreateTH1(histName, histTitle, binLabels.size(), 0, binLabels.size());
1191  // Set label names
1192  for (unsigned int i = 1; i <= binLabels.size(); i++) {
1193  histEmbeddedEventRejection->GetXaxis()->SetBinLabel(i, binLabels.at(i-1).c_str());
1194  }
1195  histEmbeddedEventRejection->GetYaxis()->SetTitle("Counts");
1196 
1197  // Rejected events in embedded event selection
1198  histName = "fHistEmbeddedEventsAttempted";
1199  histTitle = "Number of embedded events rejected by event selection before success;Number of rejected events;Counts";
1200  fHistManager.CreateTH1(histName, histTitle, 200, 0, 200);
1201 
1202  // Number of files embedded
1203  histName = "fHistNumberOfFilesEmbedded";
1204  histTitle = "Number of files which contributed events to be embedded";
1205  fHistManager.CreateTH1(histName, histTitle, 1, 0, 2);
1206 
1207  // File number which was embedded
1208  histName = "fHistAbsoluteFileNumber";
1209  histTitle = "Number of times each absolute file number was embedded";
1210  fHistManager.CreateTH1(histName, histTitle, fMaxNumberOfFiles, 0, fMaxNumberOfFiles);
1211 
1213  // Internal event cut statistics
1214  histName = "fHistInternalEventCutsStats";
1215  histTitle = "Number of events to pass each cut";
1216  binLabels = {"passedEventCuts", "centrality", "passedAllCuts"};
1217  auto histInternalEventCutsStats = fHistManager.CreateTH1(histName, histTitle, binLabels.size(), 0, binLabels.size());
1218  // Set label names
1219  for (unsigned int i = 1; i <= binLabels.size(); i++) {
1220  histInternalEventCutsStats->GetXaxis()->SetBinLabel(i, binLabels.at(i-1).c_str());
1221  }
1222  histInternalEventCutsStats->GetYaxis()->SetTitle("Number of selected events");
1223  }
1224 
1225  // Time to execute InitTree()
1226  if (fPrintTimingInfoToLog) {
1227  histName = "fInitTreeCPUtime";
1228  histTitle = "CPU time to execute InitTree() (s)";
1229  fHistManager.CreateTH1(histName, histTitle, 200, 0, 2000);
1230 
1231  histName = "fInitTreeRealtime";
1232  histTitle = "Real time to execute InitTree() (s)";
1233  fHistManager.CreateTH1(histName, histTitle, 200, 0, 2000);
1234  }
1235 
1236  // Add all histograms to output list
1237  TIter next(fHistManager.GetListOfHistograms());
1238  TObject* obj = 0;
1239  while ((obj = next())) {
1240  fOutput->Add(obj);
1241  }
1242 
1243  PostData(1, fOutput);
1244 }
1245 
1253 {
1254  // Determine which file to start with
1256 
1257  // Setup TChain
1258  fChain = new TChain(fTreeName);
1259 
1260  // Determine whether AliEn is needed
1261  for (auto filename : fFilenames)
1262  {
1263  if (filename.find("alien://") != std::string::npos) {
1264  ::ConnectToAliEn();
1265  // No point in continuing to search once we know that it is needed
1266  break;
1267  }
1268  }
1269 
1270  // Add files for TChain
1271  // See: https://stackoverflow.com/a/8533198
1272  bool wrapped = false;
1273  std::string fullPythiaXSecFilename = "";
1274  for (auto filename = fFilenames.begin() + fFilenameIndex; (filename != fFilenames.begin() + fFilenameIndex || !wrapped); filename++)
1275  {
1276  // Wraps the loop back around to the beginning
1277  if (filename == fFilenames.end())
1278  {
1279  // Explicit check is needed. Otherwise, an offset of 0 would load the 0th entry twice.
1280  if (fFilenameIndex == 0) {
1281  break;
1282  }
1283  filename = fFilenames.begin();
1284  wrapped = true;
1285  }
1286 
1287  // Add to the Chain
1288  AliDebugStream(4) << "Adding file to the embedded input chain \"" << *filename << "\".\n";
1289  fChain->Add(filename->c_str());
1290 
1291  // Determine the full pythia cross section filename based on the previously determined filename
1292  // If we have determined that it doesn't exist in the initialization then we don't repeated attempt to open
1293  // the file (which will fail)
1294  if (fPythiaXSecFilename != "") {
1295  // Could check here again whether it exists here, but almost certainly unnecessary.
1296  // Further, we won't check to ensure that rapid, repeated file access on AliEn doesn't cause any problmes!
1297  fullPythiaXSecFilename = ConstructFullPythiaXSecFilename(*filename, fPythiaXSecFilename, false);
1298 
1299  AliDebugStream(4) << "Adding pythia cross section file \"" << fullPythiaXSecFilename << "\".\n";
1300 
1301  // They will automatically be ordered the same as the files to embed!
1302  fPythiaCrossSectionFilenames.push_back(fullPythiaXSecFilename);
1303  }
1304  }
1305 
1306  // Keep track of the total number of files in the TChain to ensure that we don't start repeating within the chain
1307  fMaxNumberOfFiles = fChain->GetListOfFiles()->GetEntries();
1308 
1309  if (fFilenames.size() > fMaxNumberOfFiles) {
1310  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";
1311  }
1312 
1313  // Setup input event
1314  Bool_t res = InitEvent();
1315  if (!res) return kFALSE;
1316 
1317  return kTRUE;
1318 }
1319 
1330 std::string AliAnalysisTaskEmcalEmbeddingHelper::ConstructFullPythiaXSecFilename(std::string externalEventFilename, const std::string & pythiaFilename, bool testIfExists) const
1331 {
1332  std::string pythiaXSecFilename = "";
1333 
1334  // Remove "#*.root" if necessary
1335  if (externalEventFilename.find(".zip#") != std::string::npos) {
1336  std::size_t pos = externalEventFilename.find_last_of("#");
1337  externalEventFilename.erase(pos);
1338  }
1339 
1340  // Handle different file types
1341  if (externalEventFilename.find(".zip") != std::string::npos)
1342  {
1343  // Handle zip files
1344  pythiaXSecFilename = externalEventFilename;
1345  pythiaXSecFilename += "#";
1346  pythiaXSecFilename += pythiaFilename;
1347 
1348  // Check if the file is accessible
1349  if (testIfExists) {
1350  // Unfortunately, we cannot test for the existence of a file in an archive.
1351  // Instead, we have to tolerate TFile throwing an error (maximum of two).
1352  std::unique_ptr<TFile> fTemp(TFile::Open(pythiaXSecFilename.c_str(), "READ"));
1353 
1354  if (!fTemp) {
1355  AliDebugStream(4) << "File " << pythiaXSecFilename.c_str() << " does not exist!\n";
1356  pythiaXSecFilename = "";
1357  }
1358  else {
1359  AliDebugStream(4) << "Found pythia cross section file \"" << pythiaXSecFilename.c_str() << "\".\n";
1360  }
1361  }
1362  }
1363  else
1364  {
1365  // Handle normal root files
1366  pythiaXSecFilename = gSystem->DirName(externalEventFilename.c_str());
1367  pythiaXSecFilename += "/";
1368  pythiaXSecFilename += pythiaFilename;
1369 
1370  // Check if the file is accessible
1371  if (testIfExists) {
1372  if(::IsFileAccessible(pythiaXSecFilename)) {
1373  AliDebugStream(4) << "Found pythia cross section file \"" << pythiaXSecFilename.c_str() << "\".\n";
1374  }
1375  else {
1376  AliDebugStream(4) << "File " << pythiaXSecFilename.c_str() << " does not exist!\n";
1377  pythiaXSecFilename = "";
1378  }
1379  }
1380  }
1381 
1382  return pythiaXSecFilename;
1383 }
1384 
1391 {
1392  if (fInitializedConfiguration == false) {
1393  AliFatal("The configuration is not initialized. Check that Initialize() was called!");
1394  }
1395 
1396  // Setup TChain
1397  Bool_t res = SetupInputFiles();
1398  if (!res) { return; }
1399 
1400  // Note if getting random event access
1402  AliInfo("Random event number access enabled!");
1403  }
1404 
1405  fInitializedEmbedding = kTRUE;
1406 }
1407 
1419 {
1420  // Start the timer (for logging purposes)
1421  if (fPrintTimingInfoToLog) {
1422  fTimer.Start(kTRUE);
1423  std::cout << "InitTree() has started for file " << (fFilenameIndex + fFileNumber + 1) % fMaxNumberOfFiles << fChain->GetCurrentFile()->GetName() << "..." << std::endl;
1424  }
1425 
1426  // Load first entry of the (next) file so that we can query information about it
1427  // (it is unaccessible otherwise).
1428  // Since fUpperEntry is the total number of entries, loading it will retrieve the
1429  // next tree (in the next file) since entries are indexed starting from 0.
1430  fChain->GetEntry(fUpperEntry);
1431 
1432  // Determine tree size and current entry
1433  // Set the limits of the new tree
1435  // Fine to be += as long as we started at 0
1436  fUpperEntry += fChain->GetTree()->GetEntries();
1437 
1438  // Jump ahead at random if desired
1439  // Determines the offset into the tree
1441  TRandom3 rand(0);
1442  fOffset = TMath::Nint(rand.Rndm()*(fUpperEntry-fLowerEntry))-1;
1443  }
1444  else {
1445  fOffset = 0;
1446  }
1447 
1448  // Sets which entry to start if the try
1450 
1451  // Keep track of the number of files that we have gone through
1452  // To start from 0, we only increment if fLowerEntry > 0
1453  if (fLowerEntry > 0) {
1454  fFileNumber++;
1455  }
1456 
1457  // Add to the count the number of files which were embedded
1458  fHistManager.FillTH1("fHistNumberOfFilesEmbedded", 1);
1459  fHistManager.FillTH1("fHistAbsoluteFileNumber", (fFileNumber + fFilenameIndex) % fMaxNumberOfFiles);
1460 
1461  // Check for pythia cross section and extract if possible
1462  // fFileNumber corresponds to the next file
1463  // If there are pythia filenames, the number of match the file number of the tree.
1464  // If we previously gave up on extracting then there should be no entires
1465  if (fPythiaCrossSectionFilenames.size() > 0) {
1466  // Need to check that fFileNumber is smaller than the size of the vector because we don't check if
1469 
1470  if (!success) {
1471  AliDebugStream(3) << "Failed to retrieve cross section from xsec file. Will still attempt to get the information from the header.\n";
1472  }
1473  }
1474  else {
1475  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";
1476  }
1477  }
1478 
1479  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));
1480  // NOTE: Cannot use this print message, as it is possible that fMaxNumberOfFiles != fFilenames.size() because
1481  // invalid filenames may be included in the fFilenames count!
1482  //AliDebug(2, TString::Format("Will start embedding file %i as the %ith file beginning from entry %i.", (fFilenameIndex + fFileNumber) % fMaxNumberOfFiles, fFileNumber, fCurrentEntry));
1483 
1484  // (re)set whether we have wrapped the tree
1485  fWrappedAroundTree = false;
1486 
1487  // Note that the tree in the new file has been initialized
1488  fInitializedNewFile = kTRUE;
1489 
1490  // Stop timer (for logging purposes)
1491  if (fPrintTimingInfoToLog) {
1492  fTimer.Stop();
1493  std::cout << "InitTree() complete. CPU time: " << fTimer.CpuTime() << " (s). Real time: " << fTimer.RealTime() << " (s)." << std::endl;
1494  fHistManager.FillTH1("fInitTreeCPUtime", fTimer.CpuTime());
1495  fHistManager.FillTH1("fInitTreeRealtime", fTimer.RealTime());
1496  }
1497 
1498 }
1499 
1508 {
1509  std::unique_ptr<TFile> fxsec(TFile::Open(pythiaFileName.c_str()));
1510 
1511  if (fxsec)
1512  {
1513  int trials = 0;
1514  double crossSection = 0;
1515  double nEvents = 0;
1516  // Check if it's a tree
1517  TTree *xtree = dynamic_cast<TTree*>(fxsec->Get("Xsection"));
1518  if (xtree) {
1519  UInt_t ntrials = 0;
1520  Double_t xsection = 0;
1521  xtree->SetBranchAddress("xsection",&xsection);
1522  xtree->SetBranchAddress("ntrials",&ntrials);
1523  xtree->GetEntry(0);
1524  trials = ntrials;
1525  crossSection = xsection;
1526  // TODO: Test this on a file which has pyxsec.root!
1527  nEvents = 1.;
1528  AliFatal("Have no tested pyxsec.root files. Need to determine the proper way to get nevents!!");
1529  }
1530  else {
1531  // Check if it's instead the histograms
1532  // find the tlist we want to be independtent of the name so use the Tkey
1533  TKey* key = static_cast<TKey*>(fxsec->GetListOfKeys()->At(0));
1534  if (!key) return false;
1535  TList *list = dynamic_cast<TList*>(key->ReadObj());
1536  if (!list) return false;
1537  TProfile * crossSectionHist = static_cast<TProfile*>(list->FindObject("h1Xsec"));
1538  // check for failure
1539  if(!(crossSectionHist->GetEntries())) {
1540  // No cross seciton information available - fall back to raw
1541  AliErrorStream() << "No cross section information available in file \"" << fxsec->GetName() << "\". Will still attempt to extract cross section information from pythia header.\n";
1542  } else {
1543  // Cross section histogram filled - take it from there
1544  crossSection = crossSectionHist->GetBinContent(1);
1545  if(!crossSection) AliErrorStream() << GetName() << ": Cross section 0 for file " << pythiaFileName << std::endl;
1546  }
1547  TH1F * trialsHist = static_cast<TH1F*>(list->FindObject("h1Trials"));
1548  trials = trialsHist->GetBinContent(1);
1549  nEvents = trialsHist->GetEntries();
1550  }
1551 
1552  // If successful in retrieveing the values, normalizae the xsec and trials by the number of events
1553  // in the file. This way, we can use it as an approximate event-by-event value
1554  // We do not want to just use the overall value because some of the events may be rejected by various
1555  // event selections, so we only want that ones that were actually use. The easiest way to do so is by
1556  // filling it for each event.
1557  fPythiaTrialsFromFile = trials/nEvents;
1558  // Do __NOT__ divide by nEvents here! The value is already from a TProfile and therefore is already the mean!
1559  fPythiaCrossSectionFromFile = crossSection;
1560 
1561  return true;
1562  }
1563  else {
1564  AliDebugStream(3) << "Unable to open file \"" << pythiaFileName << "\". Will attempt to use values from the hader.";
1565  }
1566 
1567  // Could not open file
1568  return false;
1569 }
1570 
1577 {
1578  if (!fInitializedEmbedding) {
1579  AliError("Chain not initialized before running! Setting up now.");
1580  SetupEmbedding();
1581  }
1582 
1583  // Apply internal event selection
1585  fEmbeddedEventUsed = false;
1586  if (fInternalEventCuts.AcceptEvent(AliAnalysisTaskSE::InputEvent()) == true)
1587  {
1588  fEmbeddedEventUsed = true;
1589  fHistManager.FillTH1("fHistInternalEventCutsStats", "passedEventCuts", 1);
1590 
1591  // The event was accepted by AliEventCuts. Now check for additional cuts.
1592  // Centrality
1593  // NOTE: If the centrality range is the same as AliEventCuts, then simply all will pass
1594  // If a wider centrality range than in AliEventCuts is needed then it must be _entirely_
1595  // configured through manual mode.
1596  if (fCentMin != -999 && fCentMax != -999) {
1597  if (fInternalEventCuts.GetCentrality() < fCentMin || fInternalEventCuts.GetCentrality() > fCentMax) {
1598  fEmbeddedEventUsed = false;
1599  }
1600  else {
1601  fHistManager.FillTH1("fHistInternalEventCutsStats", "centrality", 1);
1602  }
1603  }
1604 
1605  if (fEmbeddedEventUsed) {
1606  // Record all cuts passed
1607  fHistManager.FillTH1("fHistInternalEventCutsStats", "passedAllCuts", 1);
1608  }
1609  }
1610 
1611  // If the internal event was rejected, then record and move on.
1612  if (fEmbeddedEventUsed == false) {
1613  if (fCreateHisto) {
1614  PostData(1, fOutput);
1615  }
1616  return;
1617  }
1618  }
1619 
1620  if (!fInitializedNewFile) {
1621  InitTree();
1622  }
1623 
1624  Bool_t res = GetNextEntry();
1625 
1626  if (!res) {
1627  AliError("Unable to get the event to embed. Nothing will be embedded.");
1628  return;
1629  }
1630 
1631  if (fCreateHisto && fOutput) {
1632  PostData(1, fOutput);
1633  }
1634 }
1635 
1640 {
1641 }
1642 
1648 {
1649  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1650  if (!mgr)
1651  {
1652  AliErrorStream() << "No analysis manager to connect to.\n";
1653  return;
1654  }
1655 
1656  // Remove the dummy task
1657  std::string dummyTaskName = GetName();
1658  dummyTaskName += "_dummyTask";
1659  TObjArray * tasks = mgr->GetTasks();
1660  if (tasks) {
1661  AliAnalysisTaskSE * dummyTask = dynamic_cast<AliAnalysisTaskSE *>(tasks->FindObject(dummyTaskName.c_str()));
1662  if (!dummyTask) {
1663  AliErrorStream() << "Could not remove dummy task \"" << dummyTaskName << "\" from analysis manager! Was it added?\n";
1664  }
1665  // Actually remove the task
1666  tasks->Remove(dummyTask);
1667  AliDebugStream(1) << "Removed dummy task named \"" << dummyTaskName << "\".\n";
1668  }
1669  else {
1670  AliErrorStream() << "Could not retrieve tasks from the analysis manager.\n";
1671  }
1672 }
1673 
1675 {
1677  return &fInternalEventCuts;
1678  }
1679  else {
1680  AliErrorStream() << "Enable manual mode in AliEventCuts (though the embedding helper) to access this object.\n";
1681  }
1682  return nullptr;
1683 }
1684 
1686 {
1687  // Get the pointer to the existing analysis manager via the static access method.
1688  //==============================================================================
1689  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1690  if (!mgr)
1691  {
1692  ::Error("AddTaskEmcalEmbeddingHelper", "No analysis manager to connect to.");
1693  return 0;
1694  }
1695 
1696  // Check the analysis type using the event handlers connected to the analysis manager.
1697  //==============================================================================
1698  AliVEventHandler* handler = mgr->GetInputEventHandler();
1699  if (!handler)
1700  {
1701  ::Error("AddTaskEmcalEmbeddingHelper", "This task requires an input event handler");
1702  return 0;
1703  }
1704 
1705  TString name = "AliAnalysisTaskEmcalEmbeddingHelper";
1706 
1707  AliAnalysisTaskEmcalEmbeddingHelper * mgrTask = static_cast<AliAnalysisTaskEmcalEmbeddingHelper *>(mgr->GetTask(name.Data()));
1708  if (mgrTask) return mgrTask;
1709 
1710  // Create the task that manages
1711  AliAnalysisTaskEmcalEmbeddingHelper * embeddingHelper = new AliAnalysisTaskEmcalEmbeddingHelper(name.Data());
1712 
1713  //-------------------------------------------------------
1714  // Final settings, pass to manager and set the containers
1715  //-------------------------------------------------------
1716 
1717  mgr->AddTask(embeddingHelper);
1718 
1719  // Create containers for input/output
1720  AliAnalysisDataContainer* cInput = mgr->GetCommonInputContainer();
1721 
1722  TString outputContainerName(name);
1723  outputContainerName += "_histos";
1724 
1725  AliAnalysisDataContainer * cOutput = mgr->CreateContainer(outputContainerName.Data(),
1726  TList::Class(),
1727  AliAnalysisManager::kOutputContainer,
1728  Form("%s", AliAnalysisManager::GetCommonFileName()));
1729 
1730  mgr->ConnectInput(embeddingHelper, 0, cInput);
1731  mgr->ConnectOutput(embeddingHelper, 1, cOutput);
1732 
1733  return embeddingHelper;
1734 }
1735 
1737 {
1738  // Get the pointer to the existing analysis manager via the static access method.
1739  //==============================================================================
1740  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1741  if (!mgr)
1742  {
1743  ::Error("ConfigureEmcalEmbeddingHelperOnLEGOTrain", "No analysis manager to connect to.");
1744  return nullptr;
1745  }
1746 
1747  // Retrieve the embedding helper
1748  auto embeddingHelperConst = AliAnalysisTaskEmcalEmbeddingHelper::GetInstance();
1749  // Cast away const-ness on the pointer since the underlying object is not const and we need to be able to modify it.
1750  auto embeddingHelper = const_cast<AliAnalysisTaskEmcalEmbeddingHelper *>(embeddingHelperConst);
1751 
1752  // Fatal if we can't find the task
1753  if (!embeddingHelper) {
1754  AliFatalClass("Could not find embedding helper, Did you remember to create it?");
1755  }
1756 
1757  AliInfoClassStream() << "Found embedding helper to configure.\n";
1758 
1759  // AliAnalysisTaskCfg will require a task to be returned, so we add a dummy task to the analysis manager,
1760  // which will be removed when the user calls Initialize(true) on the embedding helper.
1761  mgr->AddTask(new AliAnalysisTaskSE("AliAnalysisTaskEmcalEmbeddingHelper_dummyTask"));
1762 
1763  return embeddingHelper;
1764 }
1765 
1771 std::string AliAnalysisTaskEmcalEmbeddingHelper::toString(bool includeFileList) const
1772 {
1773  std::stringstream tempSS;
1774 
1775  // General embedding helper information
1776  tempSS << std::boolalpha;
1777  tempSS << GetName() << ": Embedding helper configuration:\n";
1778  tempSS << "Create histos: " << fCreateHisto << "\n";
1779  tempSS << "Pt Hard Bin: " << fPtHardBin << "\n";
1780  tempSS << "N Pt Hard Bins: " << fNPtHardBins << "\n";
1781  tempSS << "File pattern: \"" << fFilePattern << "\"\n";
1782  tempSS << "Input filename: \"" << fInputFilename << "\"\n";
1783  tempSS << "Pythia cross section filename: \"" << fPythiaXSecFilename << "\"\n";
1784  tempSS << "File list filename: \"" << fFileListFilename << "\"\n";
1785  tempSS << "Tree name: " << fTreeName << "\n";
1786  tempSS << "Print timing info to log: " << fPrintTimingInfoToLog << "\n";
1787  tempSS << "Random event number access: " << fRandomEventNumberAccess << "\n";
1788  tempSS << "Random file access: " << fRandomFileAccess << "\n";
1789  tempSS << "Starting file index: " << fFilenameIndex << "\n";
1790  tempSS << "Number of files to embed: " << fFilenames.size() << "\n";
1791  tempSS << "YAML configuration path: \"" << fConfigurationPath << "\"\n";
1792  tempSS << "Enable internal event selection: " << fUseInternalEventSelection << "\n";
1793  tempSS << "Use manual event cuts mode for internal event selection: " << fUseManualInternalEventCuts << "\n";
1794  if (fCentMin != -999 && fCentMax != -999) {
1795  tempSS << "Internal event selection centrality range: [" << fCentMin << ", " << fCentMax << "]\n";
1796  }
1797  else {
1798  tempSS << "Internal event selection centrality range disabled.\n";
1799  }
1800 
1801  std::bitset<32> triggerMask(fTriggerMask);
1802  tempSS << "\nEmbedded event settings:\n";
1803  tempSS << "Trigger mask (binary): " << triggerMask << "\n";
1804  tempSS << "Reject outliers: " << fMCRejectOutliers << "\n";
1805  tempSS << "Pt hard jet pt rejection factor: " << fPtHardJetPtRejectionFactor << "\n";
1806  tempSS << "Z vertex cut: " << fZVertexCut << "\n";
1807  tempSS << "Max difference between internal and embedded vertex: " << fMaxVertexDist << "\n";
1808 
1809  if (includeFileList) {
1810  tempSS << "\nFiles to embed:\n";
1811  for (auto filename : fFilenames) {
1812  tempSS << "\t" << filename << "\n";
1813  }
1814  }
1815 
1816  return tempSS.str();
1817 }
1818 
1826 std::ostream & AliAnalysisTaskEmcalEmbeddingHelper::Print(std::ostream & in) const {
1827  in << toString();
1828  return in;
1829 }
1830 
1839 std::ostream & operator<<(std::ostream & in, const AliAnalysisTaskEmcalEmbeddingHelper & myTask)
1840 {
1841  std::ostream & result = myTask.Print(in);
1842  return result;
1843 }
1844 
1852 {
1853  std::string temp(opt);
1854  bool includeFileList = false;
1855  if (temp == "FILELIST") {
1856  includeFileList = true;
1857  }
1858  Printf("%s", toString(includeFileList).c_str());
1859 }
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.
TStopwatch fTimer
! Timer for the InitTree() function
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.
bool fPrintTimingInfoToLog
Flag to print time to execute InitTree(), for logging purposes.
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