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