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