AliPhysics  6f1d526 (6f1d526)
 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(1),
68  fRandomEventNumberAccess(kFALSE),
69  fRandomFileAccess(kTRUE),
70  fFilePattern(""),
71  fInputFilename(""),
72  fFileListFilename(""),
73  fFilenameIndex(-1),
74  fFilenames(),
75  fTriggerMask(AliVEvent::kAny),
76  fMCRejectOutliers(false),
77  fPtHardJetPtRejectionFactor(4),
78  fZVertexCut(10),
79  fMaxVertexDist(999),
80  fExternalFile(0),
81  fCurrentEntry(0),
82  fLowerEntry(0),
83  fUpperEntry(0),
84  fOffset(0),
85  fMaxNumberOfFiles(0),
86  fFileNumber(0),
87  fInitializedConfiguration(false),
88  fInitializedEmbedding(false),
89  fInitializedNewFile(false),
90  fWrappedAroundTree(false),
91  fChain(nullptr),
92  fExternalEvent(nullptr),
93  fExternalHeader(nullptr),
94  fPythiaHeader(nullptr),
95  fPythiaTrials(0),
96  fPythiaTrialsAvg(0),
97  fPythiaCrossSection(0.),
98  fPythiaCrossSectionAvg(0.),
99  fPythiaPtHard(0.),
100  fPythiaCrossSectionFilenames(),
101  fHistManager(),
102  fOutput(nullptr)
103 {
104  if (fgInstance != 0) {
105  AliError("An instance of AliAnalysisTaskEmcalEmbeddingHelper already exists: it will be deleted!!!");
106  delete fgInstance;
107  }
108 
109  fgInstance = this;
110 }
111 
118  AliAnalysisTaskSE(name),
119  fCreateHisto(true),
120  fTreeName("aodTree"),
121  fAnchorRun(169838),
122  fPtHardBin(-1),
123  fNPtHardBins(1),
124  fRandomEventNumberAccess(kFALSE),
125  fRandomFileAccess(kTRUE),
126  fFilePattern(""),
127  fInputFilename(""),
128  fFileListFilename(""),
129  fFilenameIndex(-1),
130  fFilenames(),
131  fTriggerMask(AliVEvent::kAny),
132  fMCRejectOutliers(false),
133  fPtHardJetPtRejectionFactor(4),
134  fZVertexCut(10),
135  fMaxVertexDist(999),
136  fExternalFile(0),
137  fCurrentEntry(0),
138  fLowerEntry(0),
139  fUpperEntry(0),
140  fOffset(0),
141  fMaxNumberOfFiles(0),
142  fFileNumber(0),
143  fInitializedConfiguration(false),
144  fInitializedEmbedding(false),
145  fInitializedNewFile(false),
146  fWrappedAroundTree(false),
147  fChain(nullptr),
148  fExternalEvent(nullptr),
149  fExternalHeader(nullptr),
150  fPythiaHeader(nullptr),
151  fPythiaTrials(0),
152  fPythiaTrialsAvg(0),
153  fPythiaCrossSection(0.),
154  fPythiaCrossSectionAvg(0.),
155  fPythiaPtHard(0.),
156  fPythiaCrossSectionFilenames(),
157  fHistManager(name),
158  fOutput(nullptr)
159 {
160  if (fgInstance != 0) {
161  AliError("An instance of AliAnalysisTaskEmcalEmbeddingHelper already exists: it will be deleted!!!");
162  delete fgInstance;
163  }
164 
165  fgInstance = this;
166 
167  if (fCreateHisto) {
168  DefineOutput(1, AliEmcalList::Class());
169  }
170 }
171 
178 {
179  if (fgInstance == this) fgInstance = 0;
180  if (fExternalEvent) delete fExternalEvent;
181  if (fExternalFile) {
182  fExternalFile->Close();
183  delete fExternalFile;
184  }
185 }
186 
192 {
193  // Get file list
194  bool result = GetFilenames();
195 
196  if (result) {
198  }
199 
200  return result;
201 }
202 
227 {
228  // Determine the pattern filename if not yet set
229  if (fInputFilename == "") {
230  if (fTreeName == "aodTree") {
231  fInputFilename = "AliAOD.root";
232  }
233  else if (fTreeName == "esdTree") {
234  fInputFilename = "AliESDs.root";
235  }
236  else {
237  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()));
238  }
239  }
240 
241  // Retrieve filenames if we don't have them yet.
242  if (fFilenames.size() == 0)
243  {
244  // Handle if fPtHardBin or fAnchorRun are set
245  // This will require formatting the file pattern in the proper way to support these substitutions
246  if (fPtHardBin != -1 && fFilePattern != "") {
247  if (fAnchorRun > 0) {
248  fFilePattern = TString::Format(fFilePattern, fAnchorRun, fPtHardBin);
249  }
250  else {
251  fFilePattern = TString::Format(fFilePattern, fPtHardBin);
252  }
253  }
254 
255  // Setup AliEn access if needed
256  if (fFilePattern.Contains("alien://") || fFileListFilename.Contains("alien://")) {
257  if (!gGrid) {
258  AliInfo("Trying to connect to AliEn ...");
259  TGrid::Connect("alien://");
260  }
261  if (!gGrid) {
262  AliFatal(TString::Format("Cannot access AliEn to retrieve file list with pattern %s!", fFilePattern.Data()));
263  }
264  }
265 
266  // Retrieve AliEn filenames directly from AliEn
267  bool usedFilePattern = false;
268  if (fFilePattern.Contains("alien://")) {
269  usedFilePattern = true;
270  AliDebug(2, TString::Format("Trying to retrieve file list from AliEn with pattern file %s...", fFilePattern.Data()));
271 
272  // Create a temporary filename based on a UUID to make sure that it doesn't overwrite any files
273  if (fFileListFilename == "") {
275  }
276 
277  // The query command cannot handle "alien://" in the file pattern, so we need to remove it for the command
278  TString filePattern = fFilePattern;
279  filePattern.ReplaceAll("alien://", "");
280 
281  // Execute the grid query to get the filenames
282  AliDebug(2, TString::Format("Trying to retrieve file list from AliEn with pattern \"%s\" and input filename \"%s\"", filePattern.Data(), fInputFilename.Data()));
283  auto result = gGrid->Query(filePattern.Data(), fInputFilename.Data());
284 
285  if (result) {
286  // Loop over the result to store it in the fileList file
287  std::ofstream outFile(fFileListFilename);
288  for (int i = 0; i < result->GetEntries(); i++)
289  {
290  // "turl" corresponds to the full AliEn url
291  outFile << result->GetKey(i, "turl") << "\n";
292  }
293  outFile.close();
294  }
295  else {
296  AliErrorStream() << "Failed to run grid query\n";
297  return false;
298  }
299  }
300 
301  // Handle a filelist on AliEn
302  if (fFileListFilename.Contains("alien://")) {
303  // Check if we already used the file pattern
304  if (usedFilePattern) {
305  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";
306  }
307 
308  // Determine the local filename and copy file to local directory
309  std::string alienFilename = fFileListFilename.Data();
310  fFileListFilename = gSystem->BaseName(alienFilename.c_str());
312 
313  TFile::Cp(alienFilename.c_str(), fFileListFilename.Data());
314  }
315 
316  std::ifstream inputFile(fFileListFilename);
317 
318  // Copy available files into the filenames vector
319  // From:: https://stackoverflow.com/a/8365247
320  std::copy(std::istream_iterator<std::string>(inputFile),
321  std::istream_iterator<std::string>(),
322  std::back_inserter(fFilenames));
323 
324  inputFile.close();
325  }
326 
327  if (fFilenames.size() == 0) {
328  AliFatal(TString::Format("Filenames from pattern \"%s\" and file list \"%s\" yielded an empty list!", fFilePattern.Data(), fFileListFilename.Data()));
329  }
330 
331  // Add "#" to files in there are any zip files
332  // It is require to open the proper root file within the zip
333  for (auto filename : fFilenames)
334  {
335  if (filename.find(".zip") != std::string::npos && filename.find("#") == std::string::npos) {
336  filename += "#";
338  if (fTreeName == "aodTree") {
339  filename += "#AliAOD.root";
340  }
341  else if (fTreeName == "esdTree") {
342  filename += "#AliESDs.root";
343  }
344  else {
345  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()));
346  return false;
347  }
348  }
349  }
350 
351  AliInfoStream() << "Found " << fFilenames.size() << " files to embed\n";
352  return true;
353 }
354 
363 {
364  std::string tempStr = "";
365  if (fFileListFilename == "") {
366  tempStr = "fileList";
367  }
368  TUUID tempUUID;
369  tempStr += ".";
370  tempStr += tempUUID.AsString();
371  tempStr += ".txt";
372 
373  return tempStr;
374 }
375 
381 {
382  // This determines which file is added first to the TChain, thus determining the order of processing
383  // Random file access. Only do this if the user has no set the filename index and request random file access
384  if (fFilenameIndex == -1 && fRandomFileAccess) {
385  // - 1 ensures that we it doesn't overflow
386  fFilenameIndex = TMath::Nint(gRandom->Rndm()*fFilenames.size()) - 1;
387  // +1 to account for the fact that the filenames vector is 0 indexed.
388  AliInfo(TString::Format("Starting with random file number %i!", fFilenameIndex+1));
389  }
390  // If not random file access, then start from the beginning
391  if (fFilenameIndex >= fFilenames.size() || fFilenameIndex < 0) {
392  // Skip notifying on -1 since it will likely be set there due to constructor.
393  if (fFilenameIndex != -1) {
394  AliWarning(TString::Format("File index %i out of range from 0 to %lu! Resetting to 0!", fFilenameIndex, fFilenames.size()));
395  }
396  fFilenameIndex = 0;
397  }
398 
399  // +1 to account for the fact that the filenames vector is 0 indexed.
400  AliInfo(TString::Format("Starting with file number %i out of %lu", fFilenameIndex+1, fFilenames.size()));
401 }
402 
412 {
413  Int_t attempts = -1;
414 
415  do {
416  // Reset to start of tree
417  if (fCurrentEntry == fUpperEntry) {
419  fWrappedAroundTree = true;
420  }
421 
423  // Continue with GetEntry as normal
424  }
425  else {
426  // NOTE: On transition from one file to the next, this calls the next entry that would be expected.
427  // However, if it is for the last file, it tries to GetEntry() of one entry past the end of the last file.
428  // Normally, this would be a problem, however GetEntry() just doesn't fill the fields of an invalid index
429  // instead of throwing an error. So "invalid values" are filled for a file that doesn't exist, but then
430  // they are immediately replaced by the lines below that reset the access values and re-init the tree.
431  // The benefit of this approach is it simplies file counting (we don't need to carefully increment here
432  // and in InitTree()) and preserves the desired behavior when we are not at the last file.
433  InitTree();
434  }
435 
436  // Load current event
437  // Can be a simple less than, because fFileNumber counts from 0.
439  fChain->GetEntry(fCurrentEntry);
440  }
441  else {
442  AliError("====================================================================================================");
443  AliError("== No more files available to embed from the TChain! Restarting from the beginning of the TChain! ==");
444  AliError("== Be careful to check that this is the desired action! ==");
445  AliError("====================================================================================================");
446 
447  // Reset the relevant access values
448  // fCurrentEntry and fLowerEntry are automatically reset in InitTree()
449  fFileNumber = 0;
450  fUpperEntry = 0;
451 
452  // Re-init back to the start
453  InitTree();
454 
455  // Access the relevant entry
456  // We are certain that fFileNumber is less than fMaxNumberOfFiles, so we are resetting to start
457  fChain->GetEntry(fCurrentEntry);
458  }
459  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));
460 
461  // Set relevant event properties
463 
464  // Increment current entry
465  fCurrentEntry++;
466 
467  // Provide a check for number of attempts
468  attempts++;
469  if (attempts == 1000)
470  AliWarning("After 1000 attempts no event has been accepted by the event selection (trigger, centrality...)!");
471 
472  // Record event properties
473  if (fCreateHisto) {
475  }
476 
477  } while (!IsEventSelected());
478 
479  if (fCreateHisto) {
480  fHistManager.FillTH1("fHistEventCount", "Accepted");
481  fHistManager.FillTH1("fHistEmbeddedEventsAttempted", attempts);
482  }
483 
484  if (!fChain) return kFALSE;
485 
486  return kTRUE;
487 }
488 
494 {
495  AliDebug(4, "Set event properties");
496  fExternalHeader = fExternalEvent->GetHeader();
497 
498  // Handle pythia header if AOD
499  AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(fExternalEvent->FindListObject(AliAODMCHeader::StdBranchName()));
500  if (aodMCH) {
501  for (UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++) {
502  fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
503  if (fPythiaHeader) break;
504  }
505  }
506 
507  if (fPythiaHeader)
508  {
509  fPythiaCrossSection = fPythiaHeader->GetXsection();
510  fPythiaTrials = fPythiaHeader->Trials();
511  fPythiaPtHard = fPythiaHeader->GetPtHard();
512  // It is identically zero if the available is not available
513  if (fPythiaCrossSection == 0.) {
514  AliDebugStream(4) << "Taking the pythia cross section avg from the xsec file.\n";
516  }
517  // It is identically zero if the available is not available
518  if (fPythiaTrials == 0.) {
519  AliDebugStream(4) << "Taking the pythia trials avg from the xsec file.\n";
521  }
522  // Pt hard is inherently event-by-event and cannot by taken as a avg quantity.
523 
524  AliDebugStream(4) << "Pythia header is defined!\n";
525  AliDebugStream(4) << "fPythiaCrossSection: " << fPythiaCrossSection << "\n";
526  }
527 }
528 
533 {
534  // Fill trials, xsec, pt hard
535  fHistManager.FillTH1("fHistTrials", fPtHardBin, fPythiaTrials);
537  fHistManager.FillTH1("fHistPtHard", fPythiaPtHard);
538 }
539 
546 {
548  return kTRUE;
549  }
550 
551  if (fCreateHisto) {
552  // Keep count of number of rejected events
553  fHistManager.FillTH1("fHistEventCount", "Rejected");
554  }
555 
556  return kFALSE;
557 }
558 
565 {
566  // Physics selection
567  if (fTriggerMask != AliVEvent::kAny) {
568  UInt_t res = 0;
569  const AliESDEvent *eev = dynamic_cast<const AliESDEvent*>(InputEvent());
570  if (eev) {
571  res = (dynamic_cast<AliInputEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
572  } else {
573  const AliAODEvent *aev = dynamic_cast<const AliAODEvent*>(InputEvent());
574  if (aev) {
575  res = (dynamic_cast<AliVAODHeader*>(aev->GetHeader()))->GetOfflineTrigger();
576  }
577  }
578 
579  if ((res & fTriggerMask) == 0) {
580  AliDebug(3, Form("Event rejected due to physics selection. Event trigger mask: %d, trigger mask selection: %d.",
581  res, fTriggerMask));
582  if (fCreateHisto) {
583  fHistManager.FillTH1("fHistEmbeddedEventRejection", "PhysSel", 1);
584  }
585  return kFALSE;
586  }
587  }
588 
589  // Vertex selection
590  Double_t externalVertex[3]={0};
591  Double_t inputVertex[3]={0};
592  const AliVVertex *externalVert = fExternalEvent->GetPrimaryVertex();
593  const AliVVertex *inputVert = InputEvent()->GetPrimaryVertex();
594  if (externalVert && inputVert) {
595  externalVert->GetXYZ(externalVertex);
596  inputVert->GetXYZ(inputVertex);
597 
598  if (TMath::Abs(externalVertex[2]) > fZVertexCut) {
599  AliDebug(3, Form("Event rejected due to Z vertex selection. Event Z vertex: %f, Z vertex cut: %f",
600  externalVertex[2], fZVertexCut));
601  if (fCreateHisto) {
602  fHistManager.FillTH1("fHistEmbeddedEventRejection", "Vz", 1);
603  }
604  return kFALSE;
605  }
606  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]));
607  if (dist > fMaxVertexDist) {
608  AliDebug(3, Form("Event rejected because the distance between the current and embedded vertices is > %f. "
609  "Current event vertex (%f, %f, %f), embedded event vertex (%f, %f, %f). Distance = %f",
610  fMaxVertexDist, inputVertex[0], inputVertex[1], inputVertex[2], externalVertex[0], externalVertex[1], externalVertex[2], dist));
611  if (fCreateHisto) {
612  fHistManager.FillTH1("fHistEmbeddedEventRejection", "VertexDist", 1);
613  }
614  return kFALSE;
615  }
616  }
617 
618  // Check for pt hard bin outliers
620  {
621  // Pythia jet / pT-hard > factor
622  // This corresponds to "condition 1" in AliAnalysisTaskEmcal
623  // NOTE: The other "conditions" defined there are not really suitable to define here, since they
624  // depend on the input objects of the event
625  if (fPtHardJetPtRejectionFactor > 0.) {
626  TLorentzVector jet;
627 
628  Int_t nTriggerJets = fPythiaHeader->NTriggerJets();
629 
630  AliDebugStream(4) << "Pythia Njets: " << nTriggerJets << ", pT Hard: " << fPythiaPtHard << "\n";
631 
632  Float_t tmpjet[]={0,0,0,0};
633  for (Int_t iJet = 0; iJet< nTriggerJets; iJet++) {
634  fPythiaHeader->TriggerJet(iJet, tmpjet);
635 
636  jet.SetPxPyPzE(tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3]);
637 
638  AliDebugStream(5) << "Pythia jet " << iJet << ", pycell jet pT: " << jet.Pt() << "\n";
639 
640  //Compare jet pT and pt Hard
641  if (jet.Pt() > fPtHardJetPtRejectionFactor * fPythiaPtHard) {
642  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";
643  fHistManager.FillTH1("fHistEmbeddedEventRejection", "MCOutlier", 1);
644  return kFALSE;
645  }
646  }
647  }
648  }
649 
650  return kTRUE;
651 }
652 
659 {
660  if (!fChain) return kFALSE;
661 
662  if (!fExternalEvent) {
663  if (fTreeName == "aodTree") {
664  fExternalEvent = new AliAODEvent();
665  }
666  else if (fTreeName == "esdTree") {
667  fExternalEvent = new AliESDEvent();
668  }
669  else {
670  AliError(Form("Tree name %s not recognized!", fTreeName.Data()));
671  return kFALSE;
672  }
673  }
674 
675  fExternalEvent->ReadFromTree(fChain, fTreeName);
676 
677  return kTRUE;
678 }
679 
686 {
687  SetupEmbedding();
688 
689  if (!fCreateHisto) {
690  return;
691  }
692 
693  // Create output list
694  OpenFile(1);
695  fOutput = new AliEmcalList();
696  fOutput->SetOwner();
697 
698  // Create histograms
699  TString histName;
700  TString histTitle;
701 
702  // Cross section
703  histName = "fHistXsection";
704  histTitle = "Pythia Cross Section;p_{T} hard bin; XSection";
705  fHistManager.CreateTProfile(histName, histTitle, fNPtHardBins, 0, fNPtHardBins);
706 
707  // Trials
708  histName = "fHistTrials";
709  histTitle = "Number of Pythia Trials;p_{T} hard bin;Trials";
710  fHistManager.CreateTH1(histName, histTitle, fNPtHardBins, 0, fNPtHardBins);
711 
712  // Pt hard spectra
713  histName = "fHistPtHard";
714  histTitle = "p_{T} Hard Spectra;p_{T} hard;Counts";
715  fHistManager.CreateTH1(histName, histTitle, 500, 0, 1000);
716 
717  // Count of accepted and rejected events
718  histName = "fHistEventCount";
719  histTitle = "fHistEventCount;Result;Count";
720  auto histEventCount = fHistManager.CreateTH1(histName, histTitle, 2, 0, 2);
721  histEventCount->GetXaxis()->SetBinLabel(1,"Accepted");
722  histEventCount->GetXaxis()->SetBinLabel(2,"Rejected");
723 
724  // Event rejection reason
725  histName = "fHistEmbeddedEventRejection";
726  histTitle = "Reasons to reject embedded event";
727  std::vector<std::string> binLabels = {"PhysSel", "MCOutlier", "Vz", "VertexDist"};
728  auto fHistEmbeddedEventRejection = fHistManager.CreateTH1(histName, histTitle, binLabels.size(), 0, binLabels.size());
729  // Set label names
730  for (unsigned int i = 1; i <= binLabels.size(); i++) {
731  fHistEmbeddedEventRejection->GetXaxis()->SetBinLabel(i, binLabels.at(i-1).c_str());
732  }
733  fHistEmbeddedEventRejection->GetYaxis()->SetTitle("Counts");
734 
735  // Rejected events in embedded event selection
736  histName = "fHistEmbeddedEventsAttempted";
737  histTitle = "Number of embedded events rejected by event selection before success;Number of rejected events;Counts";
738  fHistManager.CreateTH1(histName, histTitle, 200, 0, 200);
739 
740  // Number of files embedded
741  histName = "fHistNumberOfFilesEmbedded";
742  histTitle = "Number of files which contributed events to be embeeded";
743  fHistManager.CreateTH1(histName, histTitle, 1, 0, 2);
744 
745  // Add all histograms to output list
746  TIter next(fHistManager.GetListOfHistograms());
747  TObject* obj = 0;
748  while ((obj = next())) {
749  fOutput->Add(obj);
750  }
751 
752  PostData(1, fOutput);
753 }
754 
762 {
763  // Determine which file to start with
765 
766  // Setup TChain
767  fChain = new TChain(fTreeName);
768 
769  // Determine whether AliEn is needed
770  bool requiresAlien = false;
771  for (auto filename : fFilenames)
772  {
773  if (filename.find("alien://") != std::string::npos) {
774  requiresAlien = true;
775  }
776  }
777 
778  if (requiresAlien && !gGrid) {
779  AliInfo("Trying to connect to AliEn ...");
780  TGrid::Connect("alien://");
781  }
782 
783  // Add files for TChain
784  // See: https://stackoverflow.com/a/8533198
785  bool wrapped = false;
786  TString baseFileName = "";
787  // Hanlde the pythia cross section file list
788  bool failedEntirelyToFindFile = false;
789  TString pythiaXSecFilename = "";
790  TString pythiaBaseFilename = "";
791  std::vector <std::string> pythiaBaseFilenames = {"pyxsec.root", "pyxsec_hists.root"};
792  for (auto filename = fFilenames.begin() + fFilenameIndex; (filename != fFilenames.begin() + fFilenameIndex || !wrapped); filename++)
793  {
794  // Wraps the loop back around to the beginning
795  if (filename == fFilenames.end())
796  {
797  // Explicit check is needed. Otherwise, an offset of 0 would load the 0th entry twice.
798  if (fFilenameIndex == 0) {
799  break;
800  }
801  filename = fFilenames.begin();
802  wrapped = true;
803  }
804 
805  // AccessPathName() cannot handle the "#", so we need to strip it to check that the file exists.
806  baseFileName = filename->c_str();
807  if (baseFileName.Contains(".zip#")) {
808  Ssiz_t pos = baseFileName.Last('#');
809  baseFileName.Remove(pos);
810  }
811 
812  // Ensure that the file is accessible
813  if (gSystem->AccessPathName(baseFileName)) {
814  AliError(Form("File %s does not exist! Skipping!", baseFileName.Data()));
815  // Do not process the file if it is unaccessible, but continue processing
816  continue;
817  }
818 
819  // Add to the Chain
820  AliDebugStream(4) << "Adding file to the embedded input chain \"" << filename->c_str() << "\".\n";
821  fChain->Add(filename->c_str());
822 
823  // Handle the pythia cross section (if it exists)
824  // Determiner which file it exists in (if it does exist)
825  if (pythiaBaseFilename == "" && failedEntirelyToFindFile == false) {
826  AliInfoStream() << "Attempting to determine pythia cross section filename. It can be normal to see some TFile::Init() errors!\n";
827  for (auto name : pythiaBaseFilenames) {
828  pythiaXSecFilename = DeterminePythiaXSecFilename(baseFileName, name, true);
829  if (pythiaXSecFilename != "") {
830  AliDebugStream(4) << "Found pythia cross section base filename \"" << name.c_str() << "\"\n";
831  pythiaBaseFilename = name;
832  break;
833  }
834  }
835 
836  if (pythiaBaseFilename == "") {
837  // Failed entirely - just give up on this
838  AliErrorStream() << "Failed to find pythia x sec file! Continuing with only the pythia header!\n";
839  failedEntirelyToFindFile = true;
840  }
841  else {
842  AliInfoStream() << "Found pythia cross section file \"" << pythiaBaseFilename.Data() << "\".\n";
843  }
844  }
845  // Retrieve the value based on the previously determined filename
846  // 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
847  if (failedEntirelyToFindFile == false) {
848  // Can still check whether it exists here, but we don't necessarily have to!
849  // However, we won't check to ensure that rapid file access on AliEn doesn't cause it to crash!
850  pythiaXSecFilename = DeterminePythiaXSecFilename(baseFileName, pythiaBaseFilename, false);
851 
852  AliDebugStream(4) << "Adding pythia cross section file \"" << pythiaXSecFilename.Data() << "\".\n";
853 
854  // They will automatically be ordered the same as the files to embed!
855  fPythiaCrossSectionFilenames.push_back(pythiaXSecFilename.Data());
856  }
857  }
858 
859  // Keep track of the total number of files in the TChain to ensure that we don't start repeating within the chain
860  fMaxNumberOfFiles = fChain->GetListOfFiles()->GetEntries();
861 
862  if (fFilenames.size() > fMaxNumberOfFiles) {
863  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));
864  }
865 
866  // Setup input event
867  Bool_t res = InitEvent();
868  if (!res) return kFALSE;
869 
870  return kTRUE;
871 }
872 
883 std::string AliAnalysisTaskEmcalEmbeddingHelper::DeterminePythiaXSecFilename(TString baseFileName, TString pythiaBaseFilename, bool testIfExists)
884 {
885  std::string pythiaXSecFilename = "";
886 
887  // Handle different file types
888  if (baseFileName.Contains(".zip"))
889  {
890  // Hanlde zip files
891  pythiaXSecFilename = baseFileName;
892  pythiaXSecFilename += "#";
893  pythiaXSecFilename += pythiaBaseFilename;
894 
895  // Check if the file is accessible
896  if (testIfExists) {
897  // Unfortunately, we cannot test for the existence of a file in an archive.
898  // Instead, we have to tolerate TFile throwing an error (maximum of two).
899  std::unique_ptr<TFile> fTemp(TFile::Open(pythiaXSecFilename.c_str(), "READ"));
900 
901  if (!fTemp) {
902  AliDebugStream(4) << "File " << pythiaXSecFilename.c_str() << " does not exist!\n";
903  pythiaXSecFilename = "";
904  }
905  else {
906  AliDebugStream(4) << "Found pythia cross section file \"" << pythiaXSecFilename.c_str() << "\".\n";
907  }
908  }
909  }
910  else
911  {
912  // Hanlde normal root files
913  pythiaXSecFilename = gSystem->DirName(baseFileName);
914  pythiaXSecFilename += "/";
915  pythiaXSecFilename += pythiaBaseFilename;
916 
917  // Check if the file is accessible
918  if (testIfExists) {
919  if (gSystem->AccessPathName(pythiaXSecFilename.c_str())) {
920  AliDebugStream(4) << "File " << pythiaXSecFilename.c_str() << " does not exist!\n";
921  pythiaXSecFilename = "";
922  }
923  else {
924  AliDebugStream(4) << "Found pythia cross section file \"" << pythiaXSecFilename.c_str() << "\".\n";
925  }
926  }
927  }
928 
929  return pythiaXSecFilename;
930 }
931 
938 {
939  if (fInitializedConfiguration == false) {
940  AliFatal("The configuration is not initialized. Check that Initialize() was called!");
941  }
942 
943  // Setup TChain
944  Bool_t res = SetupInputFiles();
945  if (!res) { return; }
946 
947  // Note if getting random event access
949  AliInfo("Random event number access enabled!");
950  }
951 
952  fInitializedEmbedding = kTRUE;
953 }
954 
966 {
967  // Load first entry of the (next) file so that we can query information about it
968  // (it is unaccessible otherwise).
969  // Since fUpperEntry is the total number of entries, loading it will retrieve the
970  // next tree (in the next file) since entries are indexed starting from 0.
971  fChain->GetEntry(fUpperEntry);
972 
973  // Determine tree size and current entry
974  // Set the limits of the new tree
976  // Fine to be += as long as we started at 0
977  fUpperEntry += fChain->GetTree()->GetEntries();
978 
979  // Jump ahead at random if desired
980  // Determines the offset into the tree
982  fOffset = TMath::Nint(gRandom->Rndm()*(fUpperEntry-fLowerEntry))-1;
983  }
984  else {
985  fOffset = 0;
986  }
987 
988  // Sets which entry to start if the try
990 
991  // Keep track of the number of files that we have gone through
992  // To start from 0, we only increment if fLowerEntry > 0
993  if (fLowerEntry > 0) {
994  fFileNumber++;
995  }
996 
997  // Add to the count the number of files which were embedded
998  fHistManager.FillTH1("fHistNumberOfFilesEmbedded", 1);
999 
1000  // Check for pythia cross section and extract if possible
1001  // fFileNumber corresponds to the next file
1002  // If there are pythia filenames, the number of match the file number of the tree.
1003  // If we previously gave up on extracting then there should be no entires
1004  if (fPythiaCrossSectionFilenames.size() > 0) {
1006 
1007  if (!success) {
1008  AliDebugStream(3) << "Failed to retrieve cross section from xsec file. Will still attempt to get the information from the header.\n";
1009  }
1010  }
1011 
1012  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));
1013  // NOTE: Cannot use this print message, as it is possible that fMaxNumberOfFiles != fFilenames.size() because
1014  // invalid filenames may be included in the fFilenames count!
1015  //AliDebug(2, TString::Format("Will start embedding file %i as the %ith file beginning from entry %i.", (fFilenameIndex + fFileNumber) % fMaxNumberOfFiles, fFileNumber, fCurrentEntry));
1016 
1017  // (re)set whether we have wrapped the tree
1018  fWrappedAroundTree = false;
1019 
1020  // Note that the tree in the new file has been initialized
1021  fInitializedNewFile = kTRUE;
1022 }
1023 
1032 {
1033  std::unique_ptr<TFile> fxsec(TFile::Open(pythiaFileName.c_str()));
1034 
1035  if (fxsec)
1036  {
1037  int trials = 0;
1038  double crossSection = 0;
1039  double nEvents = 0;
1040  // Check if it's a tree
1041  TTree *xtree = dynamic_cast<TTree*>(fxsec->Get("Xsection"));
1042  if (xtree) {
1043  UInt_t ntrials = 0;
1044  Double_t xsection = 0;
1045  xtree->SetBranchAddress("xsection",&xsection);
1046  xtree->SetBranchAddress("ntrials",&ntrials);
1047  xtree->GetEntry(0);
1048  trials = ntrials;
1049  crossSection = xsection;
1050  // TODO: Test this on a file which has pyxsec.root!
1051  nEvents = 1.;
1052  AliFatal("Have no tested pyxsec.root files. Need to determine the proper way to get nevents!!");
1053  }
1054  else {
1055  // Check if it's instead the histograms
1056  // find the tlist we want to be independtent of the name so use the Tkey
1057  TKey* key = static_cast<TKey*>(fxsec->GetListOfKeys()->At(0));
1058  if (!key) return false;
1059  TList *list = dynamic_cast<TList*>(key->ReadObj());
1060  if (!list) return false;
1061  TProfile * crossSectionHist = static_cast<TProfile*>(list->FindObject("h1Xsec"));
1062  // check for failure
1063  if(!(crossSectionHist->GetEntries())) {
1064  // No cross seciton information available - fall back to raw
1065  AliErrorStream() << "No cross section information available in file \"" << fxsec->GetName() << "\". Will still attempt to extract cross section information from pythia header.\n";
1066  } else {
1067  // Cross section histogram filled - take it from there
1068  crossSection = crossSectionHist->GetBinContent(1);
1069  if(!crossSection) AliErrorStream() << GetName() << ": Cross section 0 for file " << pythiaFileName << std::endl;
1070  }
1071  TH1F * trialsHist = static_cast<TH1F*>(list->FindObject("h1Trials"));
1072  trials = trialsHist->GetBinContent(1);
1073  nEvents = trialsHist->GetEntries();
1074  }
1075 
1076  // If successful in retrieveing the values, normalizae the xsec and trials by the number of events
1077  // in the file. This way, we can use it as an approximate event-by-event value
1078  // We do not want to just use the overall value because some of the events may be rejected by various
1079  // event selections, so we only want that ones that were actually use. The easiest way to do so is by
1080  // filling it for each event.
1081  fPythiaTrialsAvg = trials/nEvents;
1082  fPythiaCrossSectionAvg = crossSection/nEvents;
1083 
1084  return true;
1085  }
1086  else {
1087  AliDebugStream(3) << "Unable to open file \"" << pythiaFileName << "\". Will attempt to use values from the hader.";
1088  }
1089 
1090  // Could not open file
1091  return false;
1092 }
1093 
1100 {
1101  if (!fInitializedEmbedding) {
1102  AliError("Chain not initialized before running! Setting up now.");
1103  SetupEmbedding();
1104  }
1105 
1106  if (!fInitializedNewFile) {
1107  InitTree();
1108  }
1109 
1110  Bool_t res = GetNextEntry();
1111 
1112  if (!res) {
1113  AliError("Unable to get the event to embed. Nothing will be embedded.");
1114  return;
1115  }
1116 
1117  if (fCreateHisto && fOutput) {
1118  PostData(1, fOutput);
1119  }
1120 }
1121 
1126 {
1127 }
1128 
1137 {
1138  // Get the pointer to the existing analysis manager via the static access method.
1139  //==============================================================================
1140  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1141  if (!mgr)
1142  {
1143  ::Error("AddTaskEmcalEmbeddingHelper", "No analysis manager to connect to.");
1144  return 0;
1145  }
1146 
1147  // Check the analysis type using the event handlers connected to the analysis manager.
1148  //==============================================================================
1149  AliVEventHandler* handler = mgr->GetInputEventHandler();
1150  if (!handler)
1151  {
1152  ::Error("AddTaskEmcalEmbeddingHelper", "This task requires an input event handler");
1153  return 0;
1154  }
1155 
1156  TString name = "AliAnalysisTaskEmcalEmbeddingHelper";
1157 
1158  AliAnalysisTaskEmcalEmbeddingHelper * mgrTask = static_cast<AliAnalysisTaskEmcalEmbeddingHelper *>(mgr->GetTask(name.Data()));
1159  if (mgrTask) return mgrTask;
1160 
1161  // Create the task that manages
1162  AliAnalysisTaskEmcalEmbeddingHelper * embeddingHelper = new AliAnalysisTaskEmcalEmbeddingHelper(name.Data());
1163 
1164  //-------------------------------------------------------
1165  // Final settings, pass to manager and set the containers
1166  //-------------------------------------------------------
1167 
1168  mgr->AddTask(embeddingHelper);
1169 
1170  // Create containers for input/output
1171  AliAnalysisDataContainer* cInput = mgr->GetCommonInputContainer();
1172 
1173  TString outputContainerName(name);
1174  outputContainerName += "_histos";
1175 
1176  AliAnalysisDataContainer * cOutput = mgr->CreateContainer(outputContainerName.Data(),
1177  TList::Class(),
1178  AliAnalysisManager::kOutputContainer,
1179  Form("%s", AliAnalysisManager::GetCommonFileName()));
1180 
1181  mgr->ConnectInput(embeddingHelper, 0, cInput);
1182  mgr->ConnectOutput(embeddingHelper, 1, cOutput);
1183 
1184  return embeddingHelper;
1185 }
1186 
1192 std::string AliAnalysisTaskEmcalEmbeddingHelper::toString(bool includeFileList) const
1193 {
1194  std::stringstream tempSS;
1195 
1196  // Show the correction components
1197  tempSS << std::boolalpha;
1198  tempSS << GetName() << ": Embedding helper configuration:\n";
1199  tempSS << "Create histos: " << fCreateHisto << "\n";
1200  tempSS << "Pt Hard Bin: " << fPtHardBin << "\n";
1201  tempSS << "N Pt Hard Bins: " << fNPtHardBins << "\n";
1202  tempSS << "Anchor Run: " << fAnchorRun << "\n";
1203  tempSS << "File pattern: \"" << fFilePattern << "\"\n";
1204  tempSS << "Input filename: \"" << fInputFilename << "\"\n";
1205  tempSS << "File list filename: \"" << fFileListFilename << "\"\n";
1206  tempSS << "Tree name: " << fTreeName << "\n";
1207  tempSS << "Random event number access: " << fRandomEventNumberAccess << "\n";
1208  tempSS << "Random file access: " << fRandomFileAccess << "\n";
1209  tempSS << "Starting file index: " << fFilenameIndex << "\n";
1210  tempSS << "Number of files to embed: " << fFilenames.size() << "\n";
1211 
1212  std::bitset<32> triggerMask(fTriggerMask);
1213  tempSS << "\nEmbedded event settings:\n";
1214  tempSS << "Trigger mask (binary): " << triggerMask << "\n";
1215  tempSS << "Reject outliers: " << fMCRejectOutliers << "\n";
1216  tempSS << "Pt hard jet pt rejection factor: " << fPtHardJetPtRejectionFactor << "\n";
1217  tempSS << "Z vertex cut: " << fZVertexCut << "\n";
1218  tempSS << "Max vertex distance: " << fMaxVertexDist << "\n";
1219 
1220  if (includeFileList) {
1221  tempSS << "\nFiles to embed:\n";
1222  for (auto filename : fFilenames) {
1223  tempSS << "\t" << filename << "\n";
1224  }
1225  }
1226 
1227  return tempSS.str();
1228 }
1229 
1237 std::ostream & AliAnalysisTaskEmcalEmbeddingHelper::Print(std::ostream & in) const {
1238  in << toString();
1239  return in;
1240 }
1241 
1250 std::ostream & operator<<(std::ostream & in, const AliAnalysisTaskEmcalEmbeddingHelper & myTask)
1251 {
1252  std::ostream & result = myTask.Print(in);
1253  return result;
1254 }
1255 
1263 {
1264  std::string temp(opt);
1265  bool includeFileList = false;
1266  if (temp == "FILELIST") {
1267  includeFileList = true;
1268  }
1269  Printf("%s", toString(includeFileList).c_str());
1270 }
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
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
float Float_t
Definition: External.C:68
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_t fPtHardJetPtRejectionFactor
Factor which the pt hard bin is multiplied by to compare against pythia header jets pt...
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