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