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