AliPhysics  ec7afe5 (ec7afe5)
 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  // Keep count of the total number of events
536  fHistManager.FillTH1("fHistEventCount", "Total");
537 }
538 
545 {
546  // Trigger selection
547  if (fTriggerMask != AliVEvent::kAny) {
548  UInt_t res = 0;
549  const AliESDEvent *eev = dynamic_cast<const AliESDEvent*>(InputEvent());
550  if (eev) {
551  res = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
552  } else {
553  const AliAODEvent *aev = dynamic_cast<const AliAODEvent*>(InputEvent());
554  if (aev) {
555  res = ((AliVAODHeader*)aev->GetHeader())->GetOfflineTrigger();
556  }
557  }
558 
559  if ((res & fTriggerMask) == 0) {
560  AliDebug(3, Form("Event rejected due to physics selection. Event trigger mask: %d, trigger mask selection: %d.",
561  res, fTriggerMask));
562  return kFALSE;
563  }
564  }
565 
566  // Vertex selection
567  Double_t externalVertex[3]={0};
568  Double_t inputVertex[3]={0};
569  const AliVVertex *externalVert = fExternalEvent->GetPrimaryVertex();
570  const AliVVertex *inputVert = InputEvent()->GetPrimaryVertex();
571  if (externalVert && inputVert) {
572  externalVert->GetXYZ(externalVertex);
573  inputVert->GetXYZ(inputVertex);
574 
575  if (TMath::Abs(externalVertex[2]) > fZVertexCut) {
576  AliDebug(3, Form("Event rejected due to Z vertex selection. Event Z vertex: %f, Z vertex cut: %f",
577  externalVertex[2], fZVertexCut));
578  return kFALSE;
579  }
580  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]));
581  if (dist > fMaxVertexDist) {
582  AliDebug(3, Form("Event rejected because the distance between the current and embedded vertices is > %f. "
583  "Current event vertex (%f, %f, %f), embedded event vertex (%f, %f, %f). Distance = %f",
584  fMaxVertexDist, inputVertex[0], inputVertex[1], inputVertex[2], externalVertex[0], externalVertex[1], externalVertex[2], dist));
585  return kFALSE;
586  }
587  }
588 
589  // TODO: Can we do selection based on the contents of the external event input objects?
590  // The previous embedding task could do so by directly accessing the elements.
591  // Certainly can't do jets (say minPt of leading jet) because this has to be embedded before them.
592  // See AliJetEmbeddingFromAODTask::IsAODEventSelected()
593 
594  return kTRUE;
595 }
596 
603 {
604  if (!fChain) return kFALSE;
605 
606  if (!fExternalEvent) {
607  if (fTreeName == "aodTree") {
608  fExternalEvent = new AliAODEvent();
609  }
610  else if (fTreeName == "esdTree") {
611  fExternalEvent = new AliESDEvent();
612  }
613  else {
614  AliError(Form("Tree name %s not recognized!", fTreeName.Data()));
615  return kFALSE;
616  }
617  }
618 
619  fExternalEvent->ReadFromTree(fChain, fTreeName);
620 
621  return kTRUE;
622 }
623 
630 {
631  SetupEmbedding();
632 
633  if (!fCreateHisto) {
634  return;
635  }
636 
637  // Create output list
638  OpenFile(1);
639  fOutput = new AliEmcalList();
640  fOutput->SetOwner();
641 
642  // Create histograms
643  TString histName;
644  TString histTitle;
645 
646  // Cross section
647  histName = "fHistXsection";
648  histTitle = "Pythia Cross Section;p_{T} hard bin; XSection";
649  fHistManager.CreateTProfile(histName, histTitle, fNPtHardBins + 1, -1, fNPtHardBins);
650 
651  // Trials
652  histName = "fHistTrials";
653  histTitle = "Number of Pythia Trials;p_{T} hard bin;Trials";
654  fHistManager.CreateTH1(histName, histTitle, fNPtHardBins + 1, -1, fNPtHardBins);
655 
656  // Pt hard spectra
657  histName = "fHistPtHard";
658  histTitle = "p_{T} Hard Spectra;p_{T} hard;Counts";
659  fHistManager.CreateTH1(histName, histTitle, 500, 0, 1000);
660 
661  // Count of accepted and rejected events
662  // NOTE: This is slightly different than the one from AliAnalysisTaskEmcal due to the difficultly in
663  // properly counting the number of rejected directly. Instead, we count the total, and then
664  // rejected in just total-accepted.
665  histName = "fHistEventCount";
666  histTitle = "fHistEventCount;Result;Count";
667  auto histEventCount = fHistManager.CreateTH1(histName, histTitle, 2, 0, 2);
668  histEventCount->GetXaxis()->SetBinLabel(1,"Accepted");
669  histEventCount->GetXaxis()->SetBinLabel(2,"Total");
670 
671  // Rejected events in embedded event selection
672  histName = "fHistEmbeddingEventsRejected";
673  histTitle = "Number of embedded events rejected by event selection before success;Number of rejected events;Counts";
674  fHistManager.CreateTH1(histName, histTitle, 200, 0, 200);
675 
676  // Add all histograms to output list
677  TIter next(fHistManager.GetListOfHistograms());
678  TObject* obj = 0;
679  while ((obj = next())) {
680  fOutput->Add(obj);
681  }
682 
683  PostData(1, fOutput);
684 }
685 
693 {
694  // Determine which file to start with
696 
697  // Setup TChain
698  fChain = new TChain(fTreeName);
699 
700  // Determine whether AliEn is needed
701  bool requiresAlien = false;
702  for (auto filename : fFilenames)
703  {
704  if (filename.find("alien://") != std::string::npos) {
705  requiresAlien = true;
706  }
707  }
708 
709  if (requiresAlien && !gGrid) {
710  AliInfo("Trying to connect to AliEn ...");
711  TGrid::Connect("alien://");
712  }
713 
714  // Add files for TChain
715  // See: https://stackoverflow.com/a/8533198
716  bool wrapped = false;
717  TString baseFileName = "";
718  // Hanlde the pythia cross section file list
719  bool failedEntirelyToFindFile = false;
720  TString pythiaXSecFilename = "";
721  TString pythiaBaseFilename = "";
722  std::vector <std::string> pythiaBaseFilenames = {"pyxsec.root", "pyxsec_hists.root"};
723  for (auto filename = fFilenames.begin() + fFilenameIndex; (filename != fFilenames.begin() + fFilenameIndex || !wrapped); filename++)
724  {
725  // Wraps the loop back around to the beginning
726  if (filename == fFilenames.end())
727  {
728  // Explicit check is needed. Otherwise, an offset of 0 would load the 0th entry twice.
729  if (fFilenameIndex == 0) {
730  break;
731  }
732  filename = fFilenames.begin();
733  wrapped = true;
734  }
735 
736  // AccessPathName() cannot handle the "#", so we need to strip it to check that the file exists.
737  baseFileName = filename->c_str();
738  if (baseFileName.Contains(".zip#")) {
739  Ssiz_t pos = baseFileName.Last('#');
740  baseFileName.Remove(pos);
741  }
742 
743  // Ensure that the file is accessible
744  if (gSystem->AccessPathName(baseFileName)) {
745  AliError(Form("File %s does not exist! Skipping!", baseFileName.Data()));
746  // Do not process the file if it is unaccessible, but continue processing
747  continue;
748  }
749 
750  // Add to the Chain
751  AliDebugStream(4) << "Adding file to the embedded input chain \"" << filename->c_str() << "\".\n";
752  fChain->Add(filename->c_str());
753 
754  // Handle the pythia cross section (if it exists)
755  // Determiner which file it exists in (if it does exist)
756  if (pythiaBaseFilename == "" && failedEntirelyToFindFile == false) {
757  AliDebugStream(4) << "Attempting to determine pythia cross section filename.\n";
758  for (auto name : pythiaBaseFilenames) {
759  pythiaXSecFilename = DeterminePythiaXSecFilename(baseFileName, name, true);
760  if (pythiaXSecFilename != "") {
761  AliDebugStream(4) << "Found pythia cross section base filename \"" << name.c_str() << "\"\n";
762  pythiaBaseFilename = name;
763  break;
764  }
765  }
766 
767  if (pythiaBaseFilename == "") {
768  // Failed entirely - just give up on this
769  AliErrorStream() << "Failed to find pythia x sec file! Continuing with only the pythia header!\n";
770  failedEntirelyToFindFile = true;
771  }
772  else {
773  AliInfoStream() << "Found pythia cross section file \"" << pythiaBaseFilename.Data() << "\".\n";
774  }
775  }
776  // Retrieve the value based on the previously determined filename
777  // 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
778  if (failedEntirelyToFindFile == false) {
779  // Can still check whether it exists here, but we don't necessarily have to!
780  // However, we won't check to ensure that rapid file access on AliEn doesn't cause it to crash!
781  pythiaXSecFilename = DeterminePythiaXSecFilename(baseFileName, pythiaBaseFilename, false);
782 
783  AliDebugStream(4) << "Adding pythia cross section file \"" << pythiaXSecFilename.Data() << "\".\n";
784 
785  // They will automatically be ordered the same as the files to embed!
786  fPythiaCrossSectionFilenames.push_back(pythiaXSecFilename.Data());
787  }
788  }
789 
790  // Keep track of the total number of files in the TChain to ensure that we don't start repeating within the chain
791  fMaxNumberOfFiles = fChain->GetListOfFiles()->GetEntries();
792 
793  if (fFilenames.size() > fMaxNumberOfFiles) {
794  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));
795  }
796 
797  // Setup input event
798  Bool_t res = InitEvent();
799  if (!res) return kFALSE;
800 
801  return kTRUE;
802 }
803 
814 std::string AliAnalysisTaskEmcalEmbeddingHelper::DeterminePythiaXSecFilename(TString baseFileName, TString pythiaBaseFilename, bool testIfExists)
815 {
816  std::string pythiaXSecFilename = "";
817 
818  // Handle different file types
819  if (baseFileName.Contains(".zip"))
820  {
821  // Hanlde zip files
822  pythiaXSecFilename = baseFileName;
823  pythiaXSecFilename += "#";
824  pythiaXSecFilename += pythiaBaseFilename;
825 
826  // Check if the file is accessible
827  if (testIfExists) {
828  // Unfortunately, we cannot test for the existence of a file in an archive.
829  // Instead, we have to tolerate TFile throwing an error (maximum of two).
830  std::unique_ptr<TFile> fTemp(TFile::Open(pythiaXSecFilename.c_str(), "READ"));
831 
832  if (!fTemp) {
833  AliDebugStream(4) << "File " << pythiaXSecFilename.c_str() << " does not exist!\n";
834  pythiaXSecFilename = "";
835  }
836  else {
837  AliDebugStream(4) << "Found pythia cross section file \"" << pythiaXSecFilename.c_str() << "\".\n";
838  }
839  }
840  }
841  else
842  {
843  // Hanlde normal root files
844  pythiaXSecFilename = gSystem->DirName(baseFileName);
845  pythiaXSecFilename += "/";
846  pythiaXSecFilename += pythiaBaseFilename;
847 
848  // Check if the file is accessible
849  if (testIfExists) {
850  if (gSystem->AccessPathName(pythiaXSecFilename.c_str())) {
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 
860  return pythiaXSecFilename;
861 }
862 
869 {
870  if (fInitializedConfiguration == false) {
871  AliFatal("The configuration is not initialized. Check that Initialize() was called!");
872  }
873 
874  // Setup TChain
875  Bool_t res = SetupInputFiles();
876  if (!res) { return; }
877 
878  // Note if getting random event access
880  AliInfo("Random event number access enabled!");
881  }
882 
883  fInitializedEmbedding = kTRUE;
884 }
885 
897 {
898  // Load first entry of the (next) file so that we can query information about it
899  // (it is unaccessible otherwise).
900  // Since fUpperEntry is the total number of entries, loading it will retrieve the
901  // next tree (in the next file) since entries are indexed starting from 0.
902  fChain->GetEntry(fUpperEntry);
903 
904  // Determine tree size and current entry
905  // Set the limits of the new tree
907  // Fine to be += as long as we started at 0
908  fUpperEntry += fChain->GetTree()->GetEntries();
909 
910  // Jump ahead at random if desired
911  // Determines the offset into the tree
913  fOffset = TMath::Nint(gRandom->Rndm()*(fUpperEntry-fLowerEntry))-1;
914  }
915  else {
916  fOffset = 0;
917  }
918 
919  // Sets which entry to start if the try
921 
922  // Keep track of the number of files that we have gone through
923  // To start from 0, we only increment if fLowerEntry > 0
924  if (fLowerEntry > 0) {
925  fFileNumber++;
926  }
927 
928  // Check for pythia cross section and extract if possible
929  // fFileNumber corresponds to the next file
930  // If there are pythia filenames, the number of match the file number of the tree.
931  // If we previously gave up on extracting then there should be no entires
932  if (fPythiaCrossSectionFilenames.size() > 0) {
934 
935  if (!success) {
936  AliDebugStream(3) << "Failed to retrieve cross section from xsec file. Will still attempt to get the information from the header.\n";
937  }
938  }
939 
940  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));
941  // NOTE: Cannot use this print message, as it is possible that fMaxNumberOfFiles != fFilenames.size() because
942  // invalid filenames may be included in the fFilenames count!
943  //AliDebug(2, TString::Format("Will start embedding file %i as the %ith file beginning from entry %i.", (fFilenameIndex + fFileNumber) % fMaxNumberOfFiles, fFileNumber, fCurrentEntry));
944 
945  // (re)set whether we have wrapped the tree
946  fWrappedAroundTree = false;
947 
948  // Note that the tree in the new file has been initialized
949  fInitializedNewFile = kTRUE;
950 }
951 
960 {
961  std::unique_ptr<TFile> fxsec(TFile::Open(pythiaFileName.c_str()));
962 
963  if (fxsec)
964  {
965  int trials = 0;
966  double crossSection = 0;
967  double nEvents = 0;
968  // Check if it's a tree
969  TTree *xtree = dynamic_cast<TTree*>(fxsec->Get("Xsection"));
970  if (xtree) {
971  UInt_t ntrials = 0;
972  Double_t xsection = 0;
973  xtree->SetBranchAddress("xsection",&xsection);
974  xtree->SetBranchAddress("ntrials",&ntrials);
975  xtree->GetEntry(0);
976  trials = ntrials;
977  crossSection = xsection;
978  // TODO: Test this on a file which has pyxsec.root!
979  nEvents = 1.;
980  AliFatal("Have no tested pyxsec.root files. Need to determine the proper way to get nevents!!");
981  }
982  else {
983  // Check if it's instead the histograms
984  // find the tlist we want to be independtent of the name so use the Tkey
985  TKey* key = static_cast<TKey*>(fxsec->GetListOfKeys()->At(0));
986  if (!key) return false;
987  TList *list = dynamic_cast<TList*>(key->ReadObj());
988  if (!list) return false;
989  TProfile * crossSectionHist = static_cast<TProfile*>(list->FindObject("h1Xsec"));
990  // check for failure
991  if(!(crossSectionHist->GetEntries())) {
992  // No cross seciton information available - fall back to raw
993  AliErrorStream() << "No cross section information available in file \"" << fxsec->GetName() << "\". Will still attempt to extract cross section information from pythia header.\n";
994  } else {
995  // Cross section histogram filled - take it from there
996  crossSection = crossSectionHist->GetBinContent(1);
997  if(!crossSection) AliErrorStream() << GetName() << ": Cross section 0 for file " << pythiaFileName << std::endl;
998  }
999  TH1F * trialsHist = static_cast<TH1F*>(list->FindObject("h1Trials"));
1000  trials = trialsHist->GetBinContent(1);
1001  nEvents = trialsHist->GetEntries();
1002  }
1003 
1004  // If successful in retrieveing the values, normalizae the xsec and trials by the number of events
1005  // in the file. This way, we can use it as an approximate event-by-event value
1006  // We do not want to just use the overall value because some of the events may be rejected by various
1007  // event selections, so we only want that ones that were actually use. The easiest way to do so is by
1008  // filling it for each event.
1009  fPythiaTrialsAvg = trials/nEvents;
1010  fPythiaCrossSectionAvg = crossSection/nEvents;
1011 
1012  return true;
1013  }
1014  else {
1015  AliDebugStream(3) << "Unable to open file \"" << pythiaFileName << "\". Will attempt to use values from the hader.";
1016  }
1017 
1018  // Could not open file
1019  return false;
1020 }
1021 
1028 {
1029  if (!fInitializedEmbedding) {
1030  AliError("Chain not initialized before running! Setting up now.");
1031  SetupEmbedding();
1032  }
1033 
1034  if (!fInitializedNewFile) {
1035  InitTree();
1036  }
1037 
1038  Bool_t res = GetNextEntry();
1039 
1040  if (!res) {
1041  AliError("Unable to get the event to embed. Nothing will be embedded.");
1042  return;
1043  }
1044 
1045  if (fCreateHisto && fOutput) {
1046  PostData(1, fOutput);
1047  }
1048 }
1049 
1054 {
1055 }
1056 
1065 {
1066  // Get the pointer to the existing analysis manager via the static access method.
1067  //==============================================================================
1068  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1069  if (!mgr)
1070  {
1071  ::Error("AddTaskEmcalEmbeddingHelper", "No analysis manager to connect to.");
1072  return 0;
1073  }
1074 
1075  // Check the analysis type using the event handlers connected to the analysis manager.
1076  //==============================================================================
1077  AliVEventHandler* handler = mgr->GetInputEventHandler();
1078  if (!handler)
1079  {
1080  ::Error("AddTaskEmcalEmbeddingHelper", "This task requires an input event handler");
1081  return 0;
1082  }
1083 
1084  TString name = "AliAnalysisTaskEmcalEmbeddingHelper";
1085 
1086  AliAnalysisTaskEmcalEmbeddingHelper * mgrTask = static_cast<AliAnalysisTaskEmcalEmbeddingHelper *>(mgr->GetTask(name.Data()));
1087  if (mgrTask) return mgrTask;
1088 
1089  // Create the task that manages
1090  AliAnalysisTaskEmcalEmbeddingHelper * embeddingHelper = new AliAnalysisTaskEmcalEmbeddingHelper(name.Data());
1091 
1092  //-------------------------------------------------------
1093  // Final settings, pass to manager and set the containers
1094  //-------------------------------------------------------
1095 
1096  mgr->AddTask(embeddingHelper);
1097 
1098  // Create containers for input/output
1099  AliAnalysisDataContainer* cInput = mgr->GetCommonInputContainer();
1100 
1101  TString outputContainerName(name);
1102  outputContainerName += "_histos";
1103 
1104  AliAnalysisDataContainer * cOutput = mgr->CreateContainer(outputContainerName.Data(),
1105  TList::Class(),
1106  AliAnalysisManager::kOutputContainer,
1107  Form("%s", AliAnalysisManager::GetCommonFileName()));
1108 
1109  mgr->ConnectInput(embeddingHelper, 0, cInput);
1110  mgr->ConnectOutput(embeddingHelper, 1, cOutput);
1111 
1112  return embeddingHelper;
1113 }
1114 
1120 std::string AliAnalysisTaskEmcalEmbeddingHelper::toString(bool includeFileList) const
1121 {
1122  std::stringstream tempSS;
1123 
1124  // Show the correction components
1125  tempSS << std::boolalpha;
1126  tempSS << GetName() << ": Embedding helper configuration:\n";
1127  tempSS << "Create histos: " << fCreateHisto << "\n";
1128  tempSS << "Pt Hard Bin: " << fPtHardBin << "\n";
1129  tempSS << "N Pt Hard Bins: " << fNPtHardBins << "\n";
1130  tempSS << "Anchor Run: " << fAnchorRun << "\n";
1131  tempSS << "File pattern: \"" << fFilePattern << "\"\n";
1132  tempSS << "Input filename: \"" << fInputFilename << "\"\n";
1133  tempSS << "File list filename: \"" << fFileListFilename << "\"\n";
1134  tempSS << "Tree name: " << fTreeName << "\n";
1135  tempSS << "Random event number access: " << fRandomEventNumberAccess << "\n";
1136  tempSS << "Random file access: " << fRandomFileAccess << "\n";
1137  tempSS << "Starting file index: " << fFilenameIndex << "\n";
1138  tempSS << "Number of files to embed: " << fFilenames.size() << "\n";
1139 
1140  std::bitset<32> triggerMask(fTriggerMask);
1141  tempSS << "\nEmbedded event settings:\n";
1142  tempSS << "Trigger mask (binary): " << triggerMask << "\n";
1143  tempSS << "Z vertex cut: " << fZVertexCut << "\n";
1144  tempSS << "Max vertex distance: " << fMaxVertexDist << "\n";
1145 
1146  if (includeFileList) {
1147  tempSS << "\nFiles to embed:\n";
1148  for (auto filename : fFilenames) {
1149  tempSS << "\t" << filename << "\n";
1150  }
1151  }
1152 
1153  return tempSS.str();
1154 }
1155 
1163 std::ostream & AliAnalysisTaskEmcalEmbeddingHelper::Print(std::ostream & in) const {
1164  in << toString();
1165  return in;
1166 }
1167 
1176 std::ostream & operator<<(std::ostream & in, const AliAnalysisTaskEmcalEmbeddingHelper & myTask)
1177 {
1178  std::ostream & result = myTask.Print(in);
1179  return result;
1180 }
1181 
1189 {
1190  std::string temp(opt);
1191  bool includeFileList = false;
1192  if (temp == "FILELIST") {
1193  includeFileList = true;
1194  }
1195  Printf("%s", toString(includeFileList).c_str());
1196 }
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:658
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