AliPhysics  095eea3 (095eea3)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
37 #include <AliLog.h>
38 #include <AliAnalysisManager.h>
39 #include <AliVEvent.h>
40 #include <AliAODEvent.h>
41 #include <AliESDEvent.h>
42 #include <AliMCEvent.h>
43 #include <AliInputEventHandler.h>
44 #include <AliVHeader.h>
45 #include <AliAODMCHeader.h>
46 #include <AliGenPythiaEventHeader.h>
47 
48 #include "AliEmcalList.h"
49 
51 
55 
57 
63  fCreateHisto(true),
64  fTreeName(),
65  fAnchorRun(169838),
66  fPtHardBin(-1),
67  fNPtHardBins(0),
68  fRandomEventNumberAccess(kFALSE),
69  fRandomFileAccess(kTRUE),
70  fFilePattern(""),
71  fInputFilename(""),
72  fFileListFilename(""),
73  fFilenameIndex(-1),
74  fFilenames(),
75  fTriggerMask(AliVEvent::kAny),
76  fZVertexCut(10),
77  fMaxVertexDist(999),
78  fExternalFile(0),
79  fCurrentEntry(0),
80  fLowerEntry(0),
81  fUpperEntry(0),
82  fOffset(0),
83  fMaxNumberOfFiles(0),
84  fFileNumber(0),
85  fInitializedConfiguration(false),
86  fInitializedEmbedding(false),
87  fInitializedNewFile(false),
88  fWrappedAroundTree(false),
89  fChain(nullptr),
90  fExternalEvent(nullptr),
91  fExternalHeader(nullptr),
92  fPythiaHeader(nullptr),
93  fPythiaTrials(0),
94  fPythiaTrialsAvg(0),
95  fPythiaCrossSection(0.),
96  fPythiaCrossSectionAvg(0.),
97  fPythiaPtHard(0.),
98  fPythiaCrossSectionFilenames(),
99  fHistManager(),
100  fOutput(nullptr)
101 {
102  if (fgInstance != 0) {
103  AliError("An instance of AliAnalysisTaskEmcalEmbeddingHelper already exists: it will be deleted!!!");
104  delete fgInstance;
105  }
106 
107  fgInstance = this;
108 }
109 
116  AliAnalysisTaskSE(name),
117  fCreateHisto(true),
118  fTreeName("aodTree"),
119  fAnchorRun(169838),
120  fPtHardBin(-1),
121  fNPtHardBins(0),
122  fRandomEventNumberAccess(kFALSE),
123  fRandomFileAccess(kTRUE),
124  fFilePattern(""),
125  fInputFilename(""),
126  fFileListFilename(""),
127  fFilenameIndex(-1),
128  fFilenames(),
129  fTriggerMask(AliVEvent::kAny),
130  fZVertexCut(10),
131  fMaxVertexDist(999),
132  fExternalFile(0),
133  fCurrentEntry(0),
134  fLowerEntry(0),
135  fUpperEntry(0),
136  fOffset(0),
137  fMaxNumberOfFiles(0),
138  fFileNumber(0),
139  fInitializedConfiguration(false),
140  fInitializedEmbedding(false),
141  fInitializedNewFile(false),
142  fWrappedAroundTree(false),
143  fChain(nullptr),
144  fExternalEvent(nullptr),
145  fExternalHeader(nullptr),
146  fPythiaHeader(nullptr),
147  fPythiaTrials(0),
148  fPythiaTrialsAvg(0),
149  fPythiaCrossSection(0.),
150  fPythiaCrossSectionAvg(0.),
151  fPythiaPtHard(0.),
152  fPythiaCrossSectionFilenames(),
153  fHistManager(name),
154  fOutput(nullptr)
155 {
156  if (fgInstance != 0) {
157  AliError("An instance of AliAnalysisTaskEmcalEmbeddingHelper already exists: it will be deleted!!!");
158  delete fgInstance;
159  }
160 
161  fgInstance = this;
162 
163  if (fCreateHisto) {
164  DefineOutput(1, AliEmcalList::Class());
165  }
166 }
167 
174 {
175  if (fgInstance == this) fgInstance = 0;
176  if (fExternalEvent) delete fExternalEvent;
177  if (fExternalFile) {
178  fExternalFile->Close();
179  delete fExternalFile;
180  }
181 }
182 
188 {
189  // Get file list
190  bool result = GetFilenames();
191 
192  if (result) {
194  }
195 
196  return result;
197 }
198 
223 {
224  // Determine the pattern filename if not yet set
225  if (fInputFilename == "") {
226  if (fTreeName == "aodTree") {
227  fInputFilename = "AliAOD.root";
228  }
229  else if (fTreeName == "esdTree") {
230  fInputFilename = "AliESDs.root";
231  }
232  else {
233  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()));
234  }
235  }
236 
237  // Retrieve filenames if we don't have them yet.
238  if (fFilenames.size() == 0)
239  {
240  // Handle if fPtHardBin or fAnchorRun are set
241  // This will require formatting the file pattern in the proper way to support these substitutions
242  if (fPtHardBin != -1 && fFilePattern != "") {
243  if (fAnchorRun > 0) {
244  fFilePattern = TString::Format(fFilePattern, fAnchorRun, fPtHardBin);
245  }
246  else {
247  fFilePattern = TString::Format(fFilePattern, fPtHardBin);
248  }
249  }
250 
251  // Setup AliEn access if needed
252  if (fFilePattern.Contains("alien://") || fFileListFilename.Contains("alien://")) {
253  if (!gGrid) {
254  AliInfo("Trying to connect to AliEn ...");
255  TGrid::Connect("alien://");
256  }
257  if (!gGrid) {
258  AliFatal(TString::Format("Cannot access AliEn to retrieve file list with pattern %s!", fFilePattern.Data()));
259  }
260  }
261 
262  // Retrieve AliEn filenames directly from AliEn
263  bool usedFilePattern = false;
264  if (fFilePattern.Contains("alien://")) {
265  usedFilePattern = true;
266  AliDebug(2, TString::Format("Trying to retrieve file list from AliEn with pattern file %s...", fFilePattern.Data()));
267 
268  // Create a temporary filename based on a UUID to make sure that it doesn't overwrite any files
269  if (fFileListFilename == "") {
271  }
272 
273  // The query command cannot handle "alien://" in the file pattern, so we need to remove it for the command
274  TString filePattern = fFilePattern;
275  filePattern.ReplaceAll("alien://", "");
276 
277  // Execute the grid query to get the filenames
278  AliDebug(2, TString::Format("Trying to retrieve file list from AliEn with pattern \"%s\" and input filename \"%s\"", filePattern.Data(), fInputFilename.Data()));
279  auto result = gGrid->Query(filePattern.Data(), fInputFilename.Data());
280 
281  if (result) {
282  // Loop over the result to store it in the fileList file
283  std::ofstream outFile(fFileListFilename);
284  for (int i = 0; i < result->GetEntries(); i++)
285  {
286  // "turl" corresponds to the full AliEn url
287  outFile << result->GetKey(i, "turl") << "\n";
288  }
289  outFile.close();
290  }
291  else {
292  AliErrorStream() << "Failed to run grid query\n";
293  return false;
294  }
295  }
296 
297  // Handle a filelist on AliEn
298  if (fFileListFilename.Contains("alien://")) {
299  // Check if we already used the file pattern
300  if (usedFilePattern) {
301  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";
302  }
303 
304  // Determine the local filename and copy file to local directory
305  std::string alienFilename = fFileListFilename.Data();
306  fFileListFilename = gSystem->BaseName(alienFilename.c_str());
308 
309  TFile::Cp(alienFilename.c_str(), fFileListFilename.Data());
310  }
311 
312  std::ifstream inputFile(fFileListFilename);
313 
314  // Copy available files into the filenames vector
315  // From:: https://stackoverflow.com/a/8365247
316  std::copy(std::istream_iterator<std::string>(inputFile),
317  std::istream_iterator<std::string>(),
318  std::back_inserter(fFilenames));
319 
320  inputFile.close();
321  }
322 
323  if (fFilenames.size() == 0) {
324  AliFatal(TString::Format("Filenames from pattern \"%s\" and file list \"%s\" yielded an empty list!", fFilePattern.Data(), fFileListFilename.Data()));
325  }
326 
327  // Add "#" to files in there are any zip files
328  // It is require to open the proper root file within the zip
329  for (auto filename : fFilenames)
330  {
331  if (filename.find(".zip") != std::string::npos && filename.find("#") == std::string::npos) {
332  filename += "#";
334  if (fTreeName == "aodTree") {
335  filename += "#AliAOD.root";
336  }
337  else if (fTreeName == "esdTree") {
338  filename += "#AliESDs.root";
339  }
340  else {
341  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()));
342  return false;
343  }
344  }
345  }
346 
347  AliInfoStream() << "Found " << fFilenames.size() << " files to embed\n";
348  return true;
349 }
350 
359 {
360  std::string tempStr = "";
361  if (fFileListFilename == "") {
362  tempStr = "fileList";
363  }
364  TUUID tempUUID;
365  tempStr += ".";
366  tempStr += tempUUID.AsString();
367  tempStr += ".txt";
368 
369  return tempStr;
370 }
371 
377 {
378  // This determines which file is added first to the TChain, thus determining the order of processing
379  // Random file access. Only do this if the user has no set the filename index and request random file access
380  if (fFilenameIndex == -1 && fRandomFileAccess) {
381  // - 1 ensures that we it doesn't overflow
382  fFilenameIndex = TMath::Nint(gRandom->Rndm()*fFilenames.size()) - 1;
383  // +1 to account for the fact that the filenames vector is 0 indexed.
384  AliInfo(TString::Format("Starting with random file number %i!", fFilenameIndex+1));
385  }
386  // If not random file access, then start from the beginning
387  if (fFilenameIndex >= fFilenames.size() || fFilenameIndex < 0) {
388  // Skip notifying on -1 since it will likely be set there due to constructor.
389  if (fFilenameIndex != -1) {
390  AliWarning(TString::Format("File index %i out of range from 0 to %lu! Resetting to 0!", fFilenameIndex, fFilenames.size()));
391  }
392  fFilenameIndex = 0;
393  }
394 
395  // +1 to account for the fact that the filenames vector is 0 indexed.
396  AliInfo(TString::Format("Starting with file number %i out of %lu", fFilenameIndex+1, fFilenames.size()));
397 }
398 
408 {
409  Int_t attempts = -1;
410 
411  do {
412  // Reset to start of tree
413  if (fCurrentEntry == fUpperEntry) {
415  fWrappedAroundTree = true;
416  }
417 
419  // Continue with GetEntry as normal
420  }
421  else {
422  // NOTE: On transition from one file to the next, this calls the next entry that would be expected.
423  // However, if it is for the last file, it tries to GetEntry() of one entry past the end of the last file.
424  // Normally, this would be a problem, however GetEntry() just doesn't fill the fields of an invalid index
425  // instead of throwing an error. So "invalid values" are filled for a file that doesn't exist, but then
426  // they are immediately replaced by the lines below that reset the access values and re-init the tree.
427  // The benefit of this approach is it simplies file counting (we don't need to carefully increment here
428  // and in InitTree()) and preserves the desired behavior when we are not at the last file.
429  InitTree();
430  }
431 
432  // Load current event
433  // Can be a simple less than, because fFileNumber counts from 0.
435  fChain->GetEntry(fCurrentEntry);
436  }
437  else {
438  AliError("====================================================================================================");
439  AliError("== No more files available to embed from the TChain! Restarting from the beginning of the TChain! ==");
440  AliError("== Be careful to check that this is the desired action! ==");
441  AliError("====================================================================================================");
442 
443  // Reset the relevant access values
444  // fCurrentEntry and fLowerEntry are automatically reset in InitTree()
445  fFileNumber = 0;
446  fUpperEntry = 0;
447 
448  // Re-init back to the start
449  InitTree();
450 
451  // Access the relevant entry
452  // We are certain that fFileNumber is less than fMaxNumberOfFiles, so we are resetting to start
453  fChain->GetEntry(fCurrentEntry);
454  }
455  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));
456 
457  // Set relevant event properties
459 
460  // Increment current entry
461  fCurrentEntry++;
462 
463  // Provide a check for number of attempts
464  attempts++;
465  if (attempts == 1000)
466  AliWarning("After 1000 attempts no event has been accepted by the event selection (trigger, centrality...)!");
467 
468  // Record event properties
469  if (fCreateHisto) {
471  }
472 
473  } while (!IsEventSelected());
474 
475  if (fCreateHisto) {
476  fHistManager.FillTH1("fHistEventCount", "Accepted");
477  fHistManager.FillTH1("fHistEmbeddingEventsRejected", attempts);
478  }
479 
480  if (!fChain) return kFALSE;
481 
482  return kTRUE;
483 }
484 
490 {
491  AliDebug(4, "Set event properties");
492  fExternalHeader = fExternalEvent->GetHeader();
493 
494  // Handle pythia header if AOD
495  AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(fExternalEvent->FindListObject(AliAODMCHeader::StdBranchName()));
496  if (aodMCH) {
497  for (UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++) {
498  fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
499  if (fPythiaHeader) break;
500  }
501  }
502 
503  if (fPythiaHeader)
504  {
505  fPythiaCrossSection = fPythiaHeader->GetXsection();
506  fPythiaTrials = fPythiaHeader->Trials();
507  fPythiaPtHard = fPythiaHeader->GetPtHard();
508  // It is identically zero if the available is not available
509  if (fPythiaCrossSection == 0.) {
510  AliDebugStream(4) << "Taking the pythia cross section avg from the xsec file.\n";
512  }
513  // It is identically zero if the available is not available
514  if (fPythiaTrials == 0.) {
515  AliDebugStream(4) << "Taking the pythia trials avg from the xsec file.\n";
517  }
518  // Pt hard is inherently event-by-event and cannot by taken as a avg quantity.
519 
520  AliDebugStream(4) << "Pythia header is defined!\n";
521  AliDebugStream(4) << "fPythiaCrossSection: " << fPythiaCrossSection << "\n";
522  }
523 }
524 
529 {
530  // Fill trials, xsec, pt hard
531  fHistManager.FillTH1("fHistTrials", fPtHardBin, fPythiaTrials);
533  fHistManager.FillTH1("fHistPtHard", fPythiaPtHard);
534 }
535 
542 {
544  return kTRUE;
545  }
546 
547  if (fCreateHisto) {
548  // Keep count of number of rejected events
549  fHistManager.FillTH1("fHistEventCount", "Rejected");
550  }
551 
552  return kFALSE;
553 }
554 
561 {
562  // Trigger selection
563  if (fTriggerMask != AliVEvent::kAny) {
564  UInt_t res = 0;
565  const AliESDEvent *eev = dynamic_cast<const AliESDEvent*>(InputEvent());
566  if (eev) {
567  res = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
568  } else {
569  const AliAODEvent *aev = dynamic_cast<const AliAODEvent*>(InputEvent());
570  if (aev) {
571  res = ((AliVAODHeader*)aev->GetHeader())->GetOfflineTrigger();
572  }
573  }
574 
575  if ((res & fTriggerMask) == 0) {
576  AliDebug(3, Form("Event rejected due to physics selection. Event trigger mask: %d, trigger mask selection: %d.",
577  res, fTriggerMask));
578  return kFALSE;
579  }
580  }
581 
582  // Vertex selection
583  Double_t externalVertex[3]={0};
584  Double_t inputVertex[3]={0};
585  const AliVVertex *externalVert = fExternalEvent->GetPrimaryVertex();
586  const AliVVertex *inputVert = InputEvent()->GetPrimaryVertex();
587  if (externalVert && inputVert) {
588  externalVert->GetXYZ(externalVertex);
589  inputVert->GetXYZ(inputVertex);
590 
591  if (TMath::Abs(externalVertex[2]) > fZVertexCut) {
592  AliDebug(3, Form("Event rejected due to Z vertex selection. Event Z vertex: %f, Z vertex cut: %f",
593  externalVertex[2], fZVertexCut));
594  return kFALSE;
595  }
596  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]));
597  if (dist > fMaxVertexDist) {
598  AliDebug(3, Form("Event rejected because the distance between the current and embedded vertices is > %f. "
599  "Current event vertex (%f, %f, %f), embedded event vertex (%f, %f, %f). Distance = %f",
600  fMaxVertexDist, inputVertex[0], inputVertex[1], inputVertex[2], externalVertex[0], externalVertex[1], externalVertex[2], dist));
601  return kFALSE;
602  }
603  }
604 
605  // TODO: Can we do selection based on the contents of the external event input objects?
606  // The previous embedding task could do so by directly accessing the elements.
607  // Certainly can't do jets (say minPt of leading jet) because this has to be embedded before them.
608  // See AliJetEmbeddingFromAODTask::IsAODEventSelected()
609 
610  return kTRUE;
611 }
612 
619 {
620  if (!fChain) return kFALSE;
621 
622  if (!fExternalEvent) {
623  if (fTreeName == "aodTree") {
624  fExternalEvent = new AliAODEvent();
625  }
626  else if (fTreeName == "esdTree") {
627  fExternalEvent = new AliESDEvent();
628  }
629  else {
630  AliError(Form("Tree name %s not recognized!", fTreeName.Data()));
631  return kFALSE;
632  }
633  }
634 
635  fExternalEvent->ReadFromTree(fChain, fTreeName);
636 
637  return kTRUE;
638 }
639 
646 {
647  SetupEmbedding();
648 
649  if (!fCreateHisto) {
650  return;
651  }
652 
653  // Create output list
654  OpenFile(1);
655  fOutput = new AliEmcalList();
656  fOutput->SetOwner();
657 
658  // Create histograms
659  TString histName;
660  TString histTitle;
661 
662  // Cross section
663  histName = "fHistXsection";
664  histTitle = "Pythia Cross Section;p_{T} hard bin; XSection";
665  fHistManager.CreateTProfile(histName, histTitle, fNPtHardBins + 1, -1, fNPtHardBins);
666 
667  // Trials
668  histName = "fHistTrials";
669  histTitle = "Number of Pythia Trials;p_{T} hard bin;Trials";
670  fHistManager.CreateTH1(histName, histTitle, fNPtHardBins + 1, -1, fNPtHardBins);
671 
672  // Pt hard spectra
673  histName = "fHistPtHard";
674  histTitle = "p_{T} Hard Spectra;p_{T} hard;Counts";
675  fHistManager.CreateTH1(histName, histTitle, 500, 0, 1000);
676 
677  // Count of accepted and rejected events
678  histName = "fHistEventCount";
679  histTitle = "fHistEventCount;Result;Count";
680  auto histEventCount = fHistManager.CreateTH1(histName, histTitle, 2, 0, 2);
681  histEventCount->GetXaxis()->SetBinLabel(1,"Accepted");
682  histEventCount->GetXaxis()->SetBinLabel(2,"Rejected");
683 
684  // Rejected events in embedded event selection
685  histName = "fHistEmbeddingEventsRejected";
686  histTitle = "Number of embedded events rejected by event selection before success;Number of rejected events;Counts";
687  fHistManager.CreateTH1(histName, histTitle, 200, 0, 200);
688 
689  // Number of files embedded
690  histName = "fHistNumberOfFilesEmbedded";
691  histTitle = "Number of files which contributed events to be embeeded";
692  fHistManager.CreateTH1(histName, histTitle, 1, 0, 2);
693 
694  // Add all histograms to output list
695  TIter next(fHistManager.GetListOfHistograms());
696  TObject* obj = 0;
697  while ((obj = next())) {
698  fOutput->Add(obj);
699  }
700 
701  PostData(1, fOutput);
702 }
703 
711 {
712  // Determine which file to start with
714 
715  // Setup TChain
716  fChain = new TChain(fTreeName);
717 
718  // Determine whether AliEn is needed
719  bool requiresAlien = false;
720  for (auto filename : fFilenames)
721  {
722  if (filename.find("alien://") != std::string::npos) {
723  requiresAlien = true;
724  }
725  }
726 
727  if (requiresAlien && !gGrid) {
728  AliInfo("Trying to connect to AliEn ...");
729  TGrid::Connect("alien://");
730  }
731 
732  // Add files for TChain
733  // See: https://stackoverflow.com/a/8533198
734  bool wrapped = false;
735  TString baseFileName = "";
736  // Hanlde the pythia cross section file list
737  bool failedEntirelyToFindFile = false;
738  TString pythiaXSecFilename = "";
739  TString pythiaBaseFilename = "";
740  std::vector <std::string> pythiaBaseFilenames = {"pyxsec.root", "pyxsec_hists.root"};
741  for (auto filename = fFilenames.begin() + fFilenameIndex; (filename != fFilenames.begin() + fFilenameIndex || !wrapped); filename++)
742  {
743  // Wraps the loop back around to the beginning
744  if (filename == fFilenames.end())
745  {
746  // Explicit check is needed. Otherwise, an offset of 0 would load the 0th entry twice.
747  if (fFilenameIndex == 0) {
748  break;
749  }
750  filename = fFilenames.begin();
751  wrapped = true;
752  }
753 
754  // AccessPathName() cannot handle the "#", so we need to strip it to check that the file exists.
755  baseFileName = filename->c_str();
756  if (baseFileName.Contains(".zip#")) {
757  Ssiz_t pos = baseFileName.Last('#');
758  baseFileName.Remove(pos);
759  }
760 
761  // Ensure that the file is accessible
762  if (gSystem->AccessPathName(baseFileName)) {
763  AliError(Form("File %s does not exist! Skipping!", baseFileName.Data()));
764  // Do not process the file if it is unaccessible, but continue processing
765  continue;
766  }
767 
768  // Add to the Chain
769  AliDebugStream(4) << "Adding file to the embedded input chain \"" << filename->c_str() << "\".\n";
770  fChain->Add(filename->c_str());
771 
772  // Handle the pythia cross section (if it exists)
773  // Determiner which file it exists in (if it does exist)
774  if (pythiaBaseFilename == "" && failedEntirelyToFindFile == false) {
775  AliInfoStream() << "Attempting to determine pythia cross section filename. It can be normal to see some TFile::Init() errors!\n";
776  for (auto name : pythiaBaseFilenames) {
777  pythiaXSecFilename = DeterminePythiaXSecFilename(baseFileName, name, true);
778  if (pythiaXSecFilename != "") {
779  AliDebugStream(4) << "Found pythia cross section base filename \"" << name.c_str() << "\"\n";
780  pythiaBaseFilename = name;
781  break;
782  }
783  }
784 
785  if (pythiaBaseFilename == "") {
786  // Failed entirely - just give up on this
787  AliErrorStream() << "Failed to find pythia x sec file! Continuing with only the pythia header!\n";
788  failedEntirelyToFindFile = true;
789  }
790  else {
791  AliInfoStream() << "Found pythia cross section file \"" << pythiaBaseFilename.Data() << "\".\n";
792  }
793  }
794  // Retrieve the value based on the previously determined filename
795  // If we have determined that it doesn't exist in the first loop then we don't repeated attempt to fail to open the file
796  if (failedEntirelyToFindFile == false) {
797  // Can still check whether it exists here, but we don't necessarily have to!
798  // However, we won't check to ensure that rapid file access on AliEn doesn't cause it to crash!
799  pythiaXSecFilename = DeterminePythiaXSecFilename(baseFileName, pythiaBaseFilename, false);
800 
801  AliDebugStream(4) << "Adding pythia cross section file \"" << pythiaXSecFilename.Data() << "\".\n";
802 
803  // They will automatically be ordered the same as the files to embed!
804  fPythiaCrossSectionFilenames.push_back(pythiaXSecFilename.Data());
805  }
806  }
807 
808  // Keep track of the total number of files in the TChain to ensure that we don't start repeating within the chain
809  fMaxNumberOfFiles = fChain->GetListOfFiles()->GetEntries();
810 
811  if (fFilenames.size() > fMaxNumberOfFiles) {
812  AliWarning(TString::Format("Number of input files (%lu) is larger than the number of available files (%i). Some filenames were likely invalid!", fFilenames.size(), fMaxNumberOfFiles));
813  }
814 
815  // Setup input event
816  Bool_t res = InitEvent();
817  if (!res) return kFALSE;
818 
819  return kTRUE;
820 }
821 
832 std::string AliAnalysisTaskEmcalEmbeddingHelper::DeterminePythiaXSecFilename(TString baseFileName, TString pythiaBaseFilename, bool testIfExists)
833 {
834  std::string pythiaXSecFilename = "";
835 
836  // Handle different file types
837  if (baseFileName.Contains(".zip"))
838  {
839  // Hanlde zip files
840  pythiaXSecFilename = baseFileName;
841  pythiaXSecFilename += "#";
842  pythiaXSecFilename += pythiaBaseFilename;
843 
844  // Check if the file is accessible
845  if (testIfExists) {
846  // Unfortunately, we cannot test for the existence of a file in an archive.
847  // Instead, we have to tolerate TFile throwing an error (maximum of two).
848  std::unique_ptr<TFile> fTemp(TFile::Open(pythiaXSecFilename.c_str(), "READ"));
849 
850  if (!fTemp) {
851  AliDebugStream(4) << "File " << pythiaXSecFilename.c_str() << " does not exist!\n";
852  pythiaXSecFilename = "";
853  }
854  else {
855  AliDebugStream(4) << "Found pythia cross section file \"" << pythiaXSecFilename.c_str() << "\".\n";
856  }
857  }
858  }
859  else
860  {
861  // Hanlde normal root files
862  pythiaXSecFilename = gSystem->DirName(baseFileName);
863  pythiaXSecFilename += "/";
864  pythiaXSecFilename += pythiaBaseFilename;
865 
866  // Check if the file is accessible
867  if (testIfExists) {
868  if (gSystem->AccessPathName(pythiaXSecFilename.c_str())) {
869  AliDebugStream(4) << "File " << pythiaXSecFilename.c_str() << " does not exist!\n";
870  pythiaXSecFilename = "";
871  }
872  else {
873  AliDebugStream(4) << "Found pythia cross section file \"" << pythiaXSecFilename.c_str() << "\".\n";
874  }
875  }
876  }
877 
878  return pythiaXSecFilename;
879 }
880 
887 {
888  if (fInitializedConfiguration == false) {
889  AliFatal("The configuration is not initialized. Check that Initialize() was called!");
890  }
891 
892  // Setup TChain
893  Bool_t res = SetupInputFiles();
894  if (!res) { return; }
895 
896  // Note if getting random event access
898  AliInfo("Random event number access enabled!");
899  }
900 
901  fInitializedEmbedding = kTRUE;
902 }
903 
915 {
916  // Load first entry of the (next) file so that we can query information about it
917  // (it is unaccessible otherwise).
918  // Since fUpperEntry is the total number of entries, loading it will retrieve the
919  // next tree (in the next file) since entries are indexed starting from 0.
920  fChain->GetEntry(fUpperEntry);
921 
922  // Determine tree size and current entry
923  // Set the limits of the new tree
925  // Fine to be += as long as we started at 0
926  fUpperEntry += fChain->GetTree()->GetEntries();
927 
928  // Jump ahead at random if desired
929  // Determines the offset into the tree
931  fOffset = TMath::Nint(gRandom->Rndm()*(fUpperEntry-fLowerEntry))-1;
932  }
933  else {
934  fOffset = 0;
935  }
936 
937  // Sets which entry to start if the try
939 
940  // Keep track of the number of files that we have gone through
941  // To start from 0, we only increment if fLowerEntry > 0
942  if (fLowerEntry > 0) {
943  fFileNumber++;
944  }
945 
946  // Add to the count the number of files which were embedded
947  fHistManager.FillTH1("fHistNumberOfFilesEmbedded", 1);
948 
949  // Check for pythia cross section and extract if possible
950  // fFileNumber corresponds to the next file
951  // If there are pythia filenames, the number of match the file number of the tree.
952  // If we previously gave up on extracting then there should be no entires
953  if (fPythiaCrossSectionFilenames.size() > 0) {
955 
956  if (!success) {
957  AliDebugStream(3) << "Failed to retrieve cross section from xsec file. Will still attempt to get the information from the header.\n";
958  }
959  }
960 
961  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));
962  // NOTE: Cannot use this print message, as it is possible that fMaxNumberOfFiles != fFilenames.size() because
963  // invalid filenames may be included in the fFilenames count!
964  //AliDebug(2, TString::Format("Will start embedding file %i as the %ith file beginning from entry %i.", (fFilenameIndex + fFileNumber) % fMaxNumberOfFiles, fFileNumber, fCurrentEntry));
965 
966  // (re)set whether we have wrapped the tree
967  fWrappedAroundTree = false;
968 
969  // Note that the tree in the new file has been initialized
970  fInitializedNewFile = kTRUE;
971 }
972 
981 {
982  std::unique_ptr<TFile> fxsec(TFile::Open(pythiaFileName.c_str()));
983 
984  if (fxsec)
985  {
986  int trials = 0;
987  double crossSection = 0;
988  double nEvents = 0;
989  // Check if it's a tree
990  TTree *xtree = dynamic_cast<TTree*>(fxsec->Get("Xsection"));
991  if (xtree) {
992  UInt_t ntrials = 0;
993  Double_t xsection = 0;
994  xtree->SetBranchAddress("xsection",&xsection);
995  xtree->SetBranchAddress("ntrials",&ntrials);
996  xtree->GetEntry(0);
997  trials = ntrials;
998  crossSection = xsection;
999  // TODO: Test this on a file which has pyxsec.root!
1000  nEvents = 1.;
1001  AliFatal("Have no tested pyxsec.root files. Need to determine the proper way to get nevents!!");
1002  }
1003  else {
1004  // Check if it's instead the histograms
1005  // find the tlist we want to be independtent of the name so use the Tkey
1006  TKey* key = static_cast<TKey*>(fxsec->GetListOfKeys()->At(0));
1007  if (!key) return false;
1008  TList *list = dynamic_cast<TList*>(key->ReadObj());
1009  if (!list) return false;
1010  TProfile * crossSectionHist = static_cast<TProfile*>(list->FindObject("h1Xsec"));
1011  // check for failure
1012  if(!(crossSectionHist->GetEntries())) {
1013  // No cross seciton information available - fall back to raw
1014  AliErrorStream() << "No cross section information available in file \"" << fxsec->GetName() << "\". Will still attempt to extract cross section information from pythia header.\n";
1015  } else {
1016  // Cross section histogram filled - take it from there
1017  crossSection = crossSectionHist->GetBinContent(1);
1018  if(!crossSection) AliErrorStream() << GetName() << ": Cross section 0 for file " << pythiaFileName << std::endl;
1019  }
1020  TH1F * trialsHist = static_cast<TH1F*>(list->FindObject("h1Trials"));
1021  trials = trialsHist->GetBinContent(1);
1022  nEvents = trialsHist->GetEntries();
1023  }
1024 
1025  // If successful in retrieveing the values, normalizae the xsec and trials by the number of events
1026  // in the file. This way, we can use it as an approximate event-by-event value
1027  // We do not want to just use the overall value because some of the events may be rejected by various
1028  // event selections, so we only want that ones that were actually use. The easiest way to do so is by
1029  // filling it for each event.
1030  fPythiaTrialsAvg = trials/nEvents;
1031  fPythiaCrossSectionAvg = crossSection/nEvents;
1032 
1033  return true;
1034  }
1035  else {
1036  AliDebugStream(3) << "Unable to open file \"" << pythiaFileName << "\". Will attempt to use values from the hader.";
1037  }
1038 
1039  // Could not open file
1040  return false;
1041 }
1042 
1049 {
1050  if (!fInitializedEmbedding) {
1051  AliError("Chain not initialized before running! Setting up now.");
1052  SetupEmbedding();
1053  }
1054 
1055  if (!fInitializedNewFile) {
1056  InitTree();
1057  }
1058 
1059  Bool_t res = GetNextEntry();
1060 
1061  if (!res) {
1062  AliError("Unable to get the event to embed. Nothing will be embedded.");
1063  return;
1064  }
1065 
1066  if (fCreateHisto && fOutput) {
1067  PostData(1, fOutput);
1068  }
1069 }
1070 
1075 {
1076 }
1077 
1086 {
1087  // Get the pointer to the existing analysis manager via the static access method.
1088  //==============================================================================
1089  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1090  if (!mgr)
1091  {
1092  ::Error("AddTaskEmcalEmbeddingHelper", "No analysis manager to connect to.");
1093  return 0;
1094  }
1095 
1096  // Check the analysis type using the event handlers connected to the analysis manager.
1097  //==============================================================================
1098  AliVEventHandler* handler = mgr->GetInputEventHandler();
1099  if (!handler)
1100  {
1101  ::Error("AddTaskEmcalEmbeddingHelper", "This task requires an input event handler");
1102  return 0;
1103  }
1104 
1105  TString name = "AliAnalysisTaskEmcalEmbeddingHelper";
1106 
1107  AliAnalysisTaskEmcalEmbeddingHelper * mgrTask = static_cast<AliAnalysisTaskEmcalEmbeddingHelper *>(mgr->GetTask(name.Data()));
1108  if (mgrTask) return mgrTask;
1109 
1110  // Create the task that manages
1111  AliAnalysisTaskEmcalEmbeddingHelper * embeddingHelper = new AliAnalysisTaskEmcalEmbeddingHelper(name.Data());
1112 
1113  //-------------------------------------------------------
1114  // Final settings, pass to manager and set the containers
1115  //-------------------------------------------------------
1116 
1117  mgr->AddTask(embeddingHelper);
1118 
1119  // Create containers for input/output
1120  AliAnalysisDataContainer* cInput = mgr->GetCommonInputContainer();
1121 
1122  TString outputContainerName(name);
1123  outputContainerName += "_histos";
1124 
1125  AliAnalysisDataContainer * cOutput = mgr->CreateContainer(outputContainerName.Data(),
1126  TList::Class(),
1127  AliAnalysisManager::kOutputContainer,
1128  Form("%s", AliAnalysisManager::GetCommonFileName()));
1129 
1130  mgr->ConnectInput(embeddingHelper, 0, cInput);
1131  mgr->ConnectOutput(embeddingHelper, 1, cOutput);
1132 
1133  return embeddingHelper;
1134 }
1135 
1141 std::string AliAnalysisTaskEmcalEmbeddingHelper::toString(bool includeFileList) const
1142 {
1143  std::stringstream tempSS;
1144 
1145  // Show the correction components
1146  tempSS << std::boolalpha;
1147  tempSS << GetName() << ": Embedding helper configuration:\n";
1148  tempSS << "Create histos: " << fCreateHisto << "\n";
1149  tempSS << "Pt Hard Bin: " << fPtHardBin << "\n";
1150  tempSS << "N Pt Hard Bins: " << fNPtHardBins << "\n";
1151  tempSS << "Anchor Run: " << fAnchorRun << "\n";
1152  tempSS << "File pattern: \"" << fFilePattern << "\"\n";
1153  tempSS << "Input filename: \"" << fInputFilename << "\"\n";
1154  tempSS << "File list filename: \"" << fFileListFilename << "\"\n";
1155  tempSS << "Tree name: " << fTreeName << "\n";
1156  tempSS << "Random event number access: " << fRandomEventNumberAccess << "\n";
1157  tempSS << "Random file access: " << fRandomFileAccess << "\n";
1158  tempSS << "Starting file index: " << fFilenameIndex << "\n";
1159  tempSS << "Number of files to embed: " << fFilenames.size() << "\n";
1160 
1161  std::bitset<32> triggerMask(fTriggerMask);
1162  tempSS << "\nEmbedded event settings:\n";
1163  tempSS << "Trigger mask (binary): " << triggerMask << "\n";
1164  tempSS << "Z vertex cut: " << fZVertexCut << "\n";
1165  tempSS << "Max vertex distance: " << fMaxVertexDist << "\n";
1166 
1167  if (includeFileList) {
1168  tempSS << "\nFiles to embed:\n";
1169  for (auto filename : fFilenames) {
1170  tempSS << "\t" << filename << "\n";
1171  }
1172  }
1173 
1174  return tempSS.str();
1175 }
1176 
1184 std::ostream & AliAnalysisTaskEmcalEmbeddingHelper::Print(std::ostream & in) const {
1185  in << toString();
1186  return in;
1187 }
1188 
1197 std::ostream & operator<<(std::ostream & in, const AliAnalysisTaskEmcalEmbeddingHelper & myTask)
1198 {
1199  std::ostream & result = myTask.Print(in);
1200  return result;
1201 }
1202 
1210 {
1211  std::string temp(opt);
1212  bool includeFileList = false;
1213  if (temp == "FILELIST") {
1214  includeFileList = true;
1215  }
1216  Printf("%s", toString(includeFileList).c_str());
1217 }
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
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
Int_t fAnchorRun
Anchor run for the given pythia production.
TString fFileListFilename
Name of the file list containing paths to files to embed.
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
TFile * fExternalFile
! External file used for embedding
TList * list
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 fInitializedConfiguration
Notes if the configuration has been initialized.
AliVHeader * fExternalHeader
! Header of the current external event
AliGenPythiaEventHeader * fPythiaHeader
! Pythia header of the current external event
TRandom * gRandom
double fPythiaCrossSectionAvg
! Average pythia cross section extracted from a xsec file.
std::string DeterminePythiaXSecFilename(TString baseFileName, TString pythiaBaseFilename, bool testIfExists)
TString fInputFilename
Filename of input root files.
Int_t fMaxNumberOfFiles
! Max number of files that are in the TChain
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.
unsigned int UInt_t
Definition: External.C:33
THashList * GetListOfHistograms() const
Get the list of histograms.
Definition: THistManager.h:671
bool fInitializedEmbedding
! Notes where the TChain has been initialized for embedding
std::string toString(bool includeFileList=false) const
Implementation of task to embed external events.
TH1 * CreateTH1(const char *name, const char *title, int nbins, double xmin, double xmax, Option_t *opt="")
Create a new TH1 within the container.
int fPythiaTrialsAvg
! Average number of trials extracted from a xsec file.
void FillProfile(const char *name, double x, double y, double weight=1.)
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.
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 fPythiaPtHard
! Pt hard of the current event (extracted from the pythia header).
static AliAnalysisTaskEmcalEmbeddingHelper * AddTaskEmcalEmbeddingHelper()
Float_t nEvents[nProd]
std::ostream & operator<<(std::ostream &in, const AliAnalysisTaskEmcalEmbeddingHelper &myTask)
THistManager fHistManager
Manages access to all histograms.
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
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
const char Option_t
Definition: External.C:48
AliEmcalList * fOutput
! List which owns the output histograms to be saved
bool Bool_t
Definition: External.C:53
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.
std::vector< std::string > fFilenames
Paths to the files to embed.
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).
TList * OpenFile(const char *fname)
Definition: DrawAnaELoss.C:65
Int_t fFileNumber
! File number corresponding to the current tree