AliPhysics  ff1d528 (ff1d528)
 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  fPythiaTrialsFromFile(0),
97  fPythiaCrossSection(0.),
98  fPythiaCrossSectionFromFile(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  fPythiaTrialsFromFile(0),
153  fPythiaCrossSection(0.),
154  fPythiaCrossSectionFromFile(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) {
1005  // Need to check that fFileNumber is smaller than the size of the vector because we don't check if
1008 
1009  if (!success) {
1010  AliDebugStream(3) << "Failed to retrieve cross section from xsec file. Will still attempt to get the information from the header.\n";
1011  }
1012  }
1013  else {
1014  AliErrorStream() << "Attempted to read past the end of the pythia cross section filenames vector. File number: " << fFileNumber << ", vector size: " << fPythiaCrossSectionFilenames.size() << ".\nThis should only occur if we have run out of files to embed!\n";
1015  }
1016  }
1017 
1018  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));
1019  // NOTE: Cannot use this print message, as it is possible that fMaxNumberOfFiles != fFilenames.size() because
1020  // invalid filenames may be included in the fFilenames count!
1021  //AliDebug(2, TString::Format("Will start embedding file %i as the %ith file beginning from entry %i.", (fFilenameIndex + fFileNumber) % fMaxNumberOfFiles, fFileNumber, fCurrentEntry));
1022 
1023  // (re)set whether we have wrapped the tree
1024  fWrappedAroundTree = false;
1025 
1026  // Note that the tree in the new file has been initialized
1027  fInitializedNewFile = kTRUE;
1028 }
1029 
1038 {
1039  std::unique_ptr<TFile> fxsec(TFile::Open(pythiaFileName.c_str()));
1040 
1041  if (fxsec)
1042  {
1043  int trials = 0;
1044  double crossSection = 0;
1045  double nEvents = 0;
1046  // Check if it's a tree
1047  TTree *xtree = dynamic_cast<TTree*>(fxsec->Get("Xsection"));
1048  if (xtree) {
1049  UInt_t ntrials = 0;
1050  Double_t xsection = 0;
1051  xtree->SetBranchAddress("xsection",&xsection);
1052  xtree->SetBranchAddress("ntrials",&ntrials);
1053  xtree->GetEntry(0);
1054  trials = ntrials;
1055  crossSection = xsection;
1056  // TODO: Test this on a file which has pyxsec.root!
1057  nEvents = 1.;
1058  AliFatal("Have no tested pyxsec.root files. Need to determine the proper way to get nevents!!");
1059  }
1060  else {
1061  // Check if it's instead the histograms
1062  // find the tlist we want to be independtent of the name so use the Tkey
1063  TKey* key = static_cast<TKey*>(fxsec->GetListOfKeys()->At(0));
1064  if (!key) return false;
1065  TList *list = dynamic_cast<TList*>(key->ReadObj());
1066  if (!list) return false;
1067  TProfile * crossSectionHist = static_cast<TProfile*>(list->FindObject("h1Xsec"));
1068  // check for failure
1069  if(!(crossSectionHist->GetEntries())) {
1070  // No cross seciton information available - fall back to raw
1071  AliErrorStream() << "No cross section information available in file \"" << fxsec->GetName() << "\". Will still attempt to extract cross section information from pythia header.\n";
1072  } else {
1073  // Cross section histogram filled - take it from there
1074  crossSection = crossSectionHist->GetBinContent(1);
1075  if(!crossSection) AliErrorStream() << GetName() << ": Cross section 0 for file " << pythiaFileName << std::endl;
1076  }
1077  TH1F * trialsHist = static_cast<TH1F*>(list->FindObject("h1Trials"));
1078  trials = trialsHist->GetBinContent(1);
1079  nEvents = trialsHist->GetEntries();
1080  }
1081 
1082  // If successful in retrieveing the values, normalizae the xsec and trials by the number of events
1083  // in the file. This way, we can use it as an approximate event-by-event value
1084  // We do not want to just use the overall value because some of the events may be rejected by various
1085  // event selections, so we only want that ones that were actually use. The easiest way to do so is by
1086  // filling it for each event.
1087  fPythiaTrialsFromFile = trials/nEvents;
1088  // Do __NOT__ divide by nEvents here! The value is already from a TProfile and therefore is already the mean!
1089  fPythiaCrossSectionFromFile = crossSection;
1090 
1091  return true;
1092  }
1093  else {
1094  AliDebugStream(3) << "Unable to open file \"" << pythiaFileName << "\". Will attempt to use values from the hader.";
1095  }
1096 
1097  // Could not open file
1098  return false;
1099 }
1100 
1107 {
1108  if (!fInitializedEmbedding) {
1109  AliError("Chain not initialized before running! Setting up now.");
1110  SetupEmbedding();
1111  }
1112 
1113  if (!fInitializedNewFile) {
1114  InitTree();
1115  }
1116 
1117  Bool_t res = GetNextEntry();
1118 
1119  if (!res) {
1120  AliError("Unable to get the event to embed. Nothing will be embedded.");
1121  return;
1122  }
1123 
1124  if (fCreateHisto && fOutput) {
1125  PostData(1, fOutput);
1126  }
1127 }
1128 
1133 {
1134 }
1135 
1144 {
1145  // Get the pointer to the existing analysis manager via the static access method.
1146  //==============================================================================
1147  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1148  if (!mgr)
1149  {
1150  ::Error("AddTaskEmcalEmbeddingHelper", "No analysis manager to connect to.");
1151  return 0;
1152  }
1153 
1154  // Check the analysis type using the event handlers connected to the analysis manager.
1155  //==============================================================================
1156  AliVEventHandler* handler = mgr->GetInputEventHandler();
1157  if (!handler)
1158  {
1159  ::Error("AddTaskEmcalEmbeddingHelper", "This task requires an input event handler");
1160  return 0;
1161  }
1162 
1163  TString name = "AliAnalysisTaskEmcalEmbeddingHelper";
1164 
1165  AliAnalysisTaskEmcalEmbeddingHelper * mgrTask = static_cast<AliAnalysisTaskEmcalEmbeddingHelper *>(mgr->GetTask(name.Data()));
1166  if (mgrTask) return mgrTask;
1167 
1168  // Create the task that manages
1169  AliAnalysisTaskEmcalEmbeddingHelper * embeddingHelper = new AliAnalysisTaskEmcalEmbeddingHelper(name.Data());
1170 
1171  //-------------------------------------------------------
1172  // Final settings, pass to manager and set the containers
1173  //-------------------------------------------------------
1174 
1175  mgr->AddTask(embeddingHelper);
1176 
1177  // Create containers for input/output
1178  AliAnalysisDataContainer* cInput = mgr->GetCommonInputContainer();
1179 
1180  TString outputContainerName(name);
1181  outputContainerName += "_histos";
1182 
1183  AliAnalysisDataContainer * cOutput = mgr->CreateContainer(outputContainerName.Data(),
1184  TList::Class(),
1185  AliAnalysisManager::kOutputContainer,
1186  Form("%s", AliAnalysisManager::GetCommonFileName()));
1187 
1188  mgr->ConnectInput(embeddingHelper, 0, cInput);
1189  mgr->ConnectOutput(embeddingHelper, 1, cOutput);
1190 
1191  return embeddingHelper;
1192 }
1193 
1199 std::string AliAnalysisTaskEmcalEmbeddingHelper::toString(bool includeFileList) const
1200 {
1201  std::stringstream tempSS;
1202 
1203  // Show the correction components
1204  tempSS << std::boolalpha;
1205  tempSS << GetName() << ": Embedding helper configuration:\n";
1206  tempSS << "Create histos: " << fCreateHisto << "\n";
1207  tempSS << "Pt Hard Bin: " << fPtHardBin << "\n";
1208  tempSS << "N Pt Hard Bins: " << fNPtHardBins << "\n";
1209  tempSS << "Anchor Run: " << fAnchorRun << "\n";
1210  tempSS << "File pattern: \"" << fFilePattern << "\"\n";
1211  tempSS << "Input filename: \"" << fInputFilename << "\"\n";
1212  tempSS << "File list filename: \"" << fFileListFilename << "\"\n";
1213  tempSS << "Tree name: " << fTreeName << "\n";
1214  tempSS << "Random event number access: " << fRandomEventNumberAccess << "\n";
1215  tempSS << "Random file access: " << fRandomFileAccess << "\n";
1216  tempSS << "Starting file index: " << fFilenameIndex << "\n";
1217  tempSS << "Number of files to embed: " << fFilenames.size() << "\n";
1218 
1219  std::bitset<32> triggerMask(fTriggerMask);
1220  tempSS << "\nEmbedded event settings:\n";
1221  tempSS << "Trigger mask (binary): " << triggerMask << "\n";
1222  tempSS << "Reject outliers: " << fMCRejectOutliers << "\n";
1223  tempSS << "Pt hard jet pt rejection factor: " << fPtHardJetPtRejectionFactor << "\n";
1224  tempSS << "Z vertex cut: " << fZVertexCut << "\n";
1225  tempSS << "Max vertex distance: " << fMaxVertexDist << "\n";
1226 
1227  if (includeFileList) {
1228  tempSS << "\nFiles to embed:\n";
1229  for (auto filename : fFilenames) {
1230  tempSS << "\t" << filename << "\n";
1231  }
1232  }
1233 
1234  return tempSS.str();
1235 }
1236 
1244 std::ostream & AliAnalysisTaskEmcalEmbeddingHelper::Print(std::ostream & in) const {
1245  in << toString();
1246  return in;
1247 }
1248 
1257 std::ostream & operator<<(std::ostream & in, const AliAnalysisTaskEmcalEmbeddingHelper & myTask)
1258 {
1259  std::ostream & result = myTask.Print(in);
1260  return result;
1261 }
1262 
1270 {
1271  std::string temp(opt);
1272  bool includeFileList = false;
1273  if (temp == "FILELIST") {
1274  includeFileList = true;
1275  }
1276  Printf("%s", toString(includeFileList).c_str());
1277 }
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.
int fPythiaTrialsFromFile
! Average number of trials extracted from a xsec file.
Bool_t fRandomFileAccess
If true, it will start embedding from a random file in the input files list.
TSystem * gSystem
Int_t fCurrentEntry
! Current entry in the current tree
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
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.
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.
double fPythiaCrossSectionFromFile
! Average pythia cross section extracted from a xsec file.
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