AliPhysics  421aab4 (421aab4)
 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*>(fExternalEvent);
572  if (eev) {
573  AliFatal("Event selection is not implemented for embedding ESDs.");
574  // Unfortunately, the normal method of retrieving the trigger mask (commented out below) doesn't work for the embedded event since we don't
575  // 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
576  // probably best to avoid it if possible.
577  //
578  // Suggestions are welcome here!
579  //res = (dynamic_cast<AliInputEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
580  } else {
581  const AliAODEvent *aev = dynamic_cast<const AliAODEvent*>(fExternalEvent);
582  if (aev) {
583  res = (dynamic_cast<AliVAODHeader*>(aev->GetHeader()))->GetOfflineTrigger();
584  }
585  }
586 
587  if ((res & fTriggerMask) == 0) {
588  AliDebug(3, Form("Event rejected due to physics selection. Event trigger mask: %d, trigger mask selection: %d.",
589  res, fTriggerMask));
590  if (fCreateHisto) {
591  fHistManager.FillTH1("fHistEmbeddedEventRejection", "PhysSel", 1);
592  }
593  return kFALSE;
594  }
595  }
596 
597  // Vertex selection
598  Double_t externalVertex[3]={0};
599  Double_t inputVertex[3]={0};
600  const AliVVertex *externalVert = fExternalEvent->GetPrimaryVertex();
601  const AliVVertex *inputVert = AliAnalysisTaskSE::InputEvent()->GetPrimaryVertex();
602  if (externalVert && inputVert) {
603  externalVert->GetXYZ(externalVertex);
604  inputVert->GetXYZ(inputVertex);
605 
606  if (TMath::Abs(externalVertex[2]) > fZVertexCut) {
607  AliDebug(3, Form("Event rejected due to Z vertex selection. Event Z vertex: %f, Z vertex cut: %f",
608  externalVertex[2], fZVertexCut));
609  if (fCreateHisto) {
610  fHistManager.FillTH1("fHistEmbeddedEventRejection", "Vz", 1);
611  }
612  return kFALSE;
613  }
614  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]));
615  if (dist > fMaxVertexDist) {
616  AliDebug(3, Form("Event rejected because the distance between the current and embedded vertices is > %f. "
617  "Current event vertex (%f, %f, %f), embedded event vertex (%f, %f, %f). Distance = %f",
618  fMaxVertexDist, inputVertex[0], inputVertex[1], inputVertex[2], externalVertex[0], externalVertex[1], externalVertex[2], dist));
619  if (fCreateHisto) {
620  fHistManager.FillTH1("fHistEmbeddedEventRejection", "VertexDist", 1);
621  }
622  return kFALSE;
623  }
624  }
625 
626  // Check for pt hard bin outliers
628  {
629  // Pythia jet / pT-hard > factor
630  // This corresponds to "condition 1" in AliAnalysisTaskEmcal
631  // NOTE: The other "conditions" defined there are not really suitable to define here, since they
632  // depend on the input objects of the event
633  if (fPtHardJetPtRejectionFactor > 0.) {
634  TLorentzVector jet;
635 
636  Int_t nTriggerJets = fPythiaHeader->NTriggerJets();
637 
638  AliDebugStream(4) << "Pythia Njets: " << nTriggerJets << ", pT Hard: " << fPythiaPtHard << "\n";
639 
640  Float_t tmpjet[]={0,0,0,0};
641  for (Int_t iJet = 0; iJet< nTriggerJets; iJet++) {
642  fPythiaHeader->TriggerJet(iJet, tmpjet);
643 
644  jet.SetPxPyPzE(tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3]);
645 
646  AliDebugStream(5) << "Pythia jet " << iJet << ", pycell jet pT: " << jet.Pt() << "\n";
647 
648  //Compare jet pT and pt Hard
649  if (jet.Pt() > fPtHardJetPtRejectionFactor * fPythiaPtHard) {
650  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";
651  fHistManager.FillTH1("fHistEmbeddedEventRejection", "MCOutlier", 1);
652  return kFALSE;
653  }
654  }
655  }
656  }
657 
658  return kTRUE;
659 }
660 
667 {
668  if (!fChain) return kFALSE;
669 
670  if (!fExternalEvent) {
671  if (fTreeName == "aodTree") {
672  fExternalEvent = new AliAODEvent();
673  }
674  else if (fTreeName == "esdTree") {
675  fExternalEvent = new AliESDEvent();
676  }
677  else {
678  AliError(Form("Tree name %s not recognized!", fTreeName.Data()));
679  return kFALSE;
680  }
681  }
682 
683  fExternalEvent->ReadFromTree(fChain, fTreeName);
684 
685  return kTRUE;
686 }
687 
694 {
695  SetupEmbedding();
696 
697  if (!fCreateHisto) {
698  return;
699  }
700 
701  // Create output list
702  OpenFile(1);
703  fOutput = new AliEmcalList();
704  fOutput->SetOwner();
705 
706  // Create histograms
707  TString histName;
708  TString histTitle;
709 
710  // Cross section
711  histName = "fHistXsection";
712  histTitle = "Pythia Cross Section;p_{T} hard bin; XSection";
713  fHistManager.CreateTProfile(histName, histTitle, fNPtHardBins, 0, fNPtHardBins);
714 
715  // Trials
716  histName = "fHistTrials";
717  histTitle = "Number of Pythia Trials;p_{T} hard bin;Trials";
718  fHistManager.CreateTH1(histName, histTitle, fNPtHardBins, 0, fNPtHardBins);
719 
720  // Pt hard spectra
721  histName = "fHistPtHard";
722  histTitle = "p_{T} Hard Spectra;p_{T} hard;Counts";
723  fHistManager.CreateTH1(histName, histTitle, 500, 0, 1000);
724 
725  // Count of accepted and rejected events
726  histName = "fHistEventCount";
727  histTitle = "fHistEventCount;Result;Count";
728  auto histEventCount = fHistManager.CreateTH1(histName, histTitle, 2, 0, 2);
729  histEventCount->GetXaxis()->SetBinLabel(1,"Accepted");
730  histEventCount->GetXaxis()->SetBinLabel(2,"Rejected");
731 
732  // Event rejection reason
733  histName = "fHistEmbeddedEventRejection";
734  histTitle = "Reasons to reject embedded event";
735  std::vector<std::string> binLabels = {"PhysSel", "MCOutlier", "Vz", "VertexDist"};
736  auto fHistEmbeddedEventRejection = fHistManager.CreateTH1(histName, histTitle, binLabels.size(), 0, binLabels.size());
737  // Set label names
738  for (unsigned int i = 1; i <= binLabels.size(); i++) {
739  fHistEmbeddedEventRejection->GetXaxis()->SetBinLabel(i, binLabels.at(i-1).c_str());
740  }
741  fHistEmbeddedEventRejection->GetYaxis()->SetTitle("Counts");
742 
743  // Rejected events in embedded event selection
744  histName = "fHistEmbeddedEventsAttempted";
745  histTitle = "Number of embedded events rejected by event selection before success;Number of rejected events;Counts";
746  fHistManager.CreateTH1(histName, histTitle, 200, 0, 200);
747 
748  // Number of files embedded
749  histName = "fHistNumberOfFilesEmbedded";
750  histTitle = "Number of files which contributed events to be embedded";
751  fHistManager.CreateTH1(histName, histTitle, 1, 0, 2);
752 
753  // File number which was embedded
754  histName = "fHistAbsoluteFileNumber";
755  histTitle = "Number of times each absolute file number was embedded";
756  fHistManager.CreateTH1(histName, histTitle, fMaxNumberOfFiles, 0, fMaxNumberOfFiles);
757 
758  // Add all histograms to output list
759  TIter next(fHistManager.GetListOfHistograms());
760  TObject* obj = 0;
761  while ((obj = next())) {
762  fOutput->Add(obj);
763  }
764 
765  PostData(1, fOutput);
766 }
767 
775 {
776  // Determine which file to start with
778 
779  // Setup TChain
780  fChain = new TChain(fTreeName);
781 
782  // Determine whether AliEn is needed
783  bool requiresAlien = false;
784  for (auto filename : fFilenames)
785  {
786  if (filename.find("alien://") != std::string::npos) {
787  requiresAlien = true;
788  }
789  }
790 
791  if (requiresAlien && !gGrid) {
792  AliInfo("Trying to connect to AliEn ...");
793  TGrid::Connect("alien://");
794  }
795 
796  // Add files for TChain
797  // See: https://stackoverflow.com/a/8533198
798  bool wrapped = false;
799  TString baseFileName = "";
800  // Hanlde the pythia cross section file list
801  bool failedEntirelyToFindFile = false;
802  TString pythiaXSecFilename = "";
803  TString pythiaBaseFilename = "";
804  std::vector <std::string> pythiaBaseFilenames = {"pyxsec.root", "pyxsec_hists.root"};
805  for (auto filename = fFilenames.begin() + fFilenameIndex; (filename != fFilenames.begin() + fFilenameIndex || !wrapped); filename++)
806  {
807  // Wraps the loop back around to the beginning
808  if (filename == fFilenames.end())
809  {
810  // Explicit check is needed. Otherwise, an offset of 0 would load the 0th entry twice.
811  if (fFilenameIndex == 0) {
812  break;
813  }
814  filename = fFilenames.begin();
815  wrapped = true;
816  }
817 
818  // AccessPathName() cannot handle the "#", so we need to strip it to check that the file exists.
819  baseFileName = filename->c_str();
820  if (baseFileName.Contains(".zip#")) {
821  Ssiz_t pos = baseFileName.Last('#');
822  baseFileName.Remove(pos);
823  }
824 
825  // Ensure that the file is accessible
826  if (gSystem->AccessPathName(baseFileName)) {
827  AliError(Form("File %s does not exist! Skipping!", baseFileName.Data()));
828  // Do not process the file if it is unaccessible, but continue processing
829  continue;
830  }
831 
832  // Add to the Chain
833  AliDebugStream(4) << "Adding file to the embedded input chain \"" << filename->c_str() << "\".\n";
834  fChain->Add(filename->c_str());
835 
836  // Handle the pythia cross section (if it exists)
837  // Determiner which file it exists in (if it does exist)
838  if (pythiaBaseFilename == "" && failedEntirelyToFindFile == false) {
839  AliInfoStream() << "Attempting to determine pythia cross section filename. It can be normal to see some TFile::Init() errors!\n";
840  for (auto name : pythiaBaseFilenames) {
841  pythiaXSecFilename = DeterminePythiaXSecFilename(baseFileName, name, true);
842  if (pythiaXSecFilename != "") {
843  AliDebugStream(4) << "Found pythia cross section base filename \"" << name.c_str() << "\"\n";
844  pythiaBaseFilename = name;
845  break;
846  }
847  }
848 
849  if (pythiaBaseFilename == "") {
850  // Failed entirely - just give up on this
851  AliErrorStream() << "Failed to find pythia x sec file! Continuing with only the pythia header!\n";
852  failedEntirelyToFindFile = true;
853  }
854  else {
855  AliInfoStream() << "Found pythia cross section file \"" << pythiaBaseFilename.Data() << "\".\n";
856  }
857  }
858  // Retrieve the value based on the previously determined filename
859  // 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
860  if (failedEntirelyToFindFile == false) {
861  // Can still check whether it exists here, but we don't necessarily have to!
862  // However, we won't check to ensure that rapid file access on AliEn doesn't cause it to crash!
863  pythiaXSecFilename = DeterminePythiaXSecFilename(baseFileName, pythiaBaseFilename, false);
864 
865  AliDebugStream(4) << "Adding pythia cross section file \"" << pythiaXSecFilename.Data() << "\".\n";
866 
867  // They will automatically be ordered the same as the files to embed!
868  fPythiaCrossSectionFilenames.push_back(pythiaXSecFilename.Data());
869  }
870  }
871 
872  // Keep track of the total number of files in the TChain to ensure that we don't start repeating within the chain
873  fMaxNumberOfFiles = fChain->GetListOfFiles()->GetEntries();
874 
875  if (fFilenames.size() > fMaxNumberOfFiles) {
876  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));
877  }
878 
879  // Setup input event
880  Bool_t res = InitEvent();
881  if (!res) return kFALSE;
882 
883  return kTRUE;
884 }
885 
896 std::string AliAnalysisTaskEmcalEmbeddingHelper::DeterminePythiaXSecFilename(TString baseFileName, TString pythiaBaseFilename, bool testIfExists)
897 {
898  std::string pythiaXSecFilename = "";
899 
900  // Handle different file types
901  if (baseFileName.Contains(".zip"))
902  {
903  // Hanlde zip files
904  pythiaXSecFilename = baseFileName;
905  pythiaXSecFilename += "#";
906  pythiaXSecFilename += pythiaBaseFilename;
907 
908  // Check if the file is accessible
909  if (testIfExists) {
910  // Unfortunately, we cannot test for the existence of a file in an archive.
911  // Instead, we have to tolerate TFile throwing an error (maximum of two).
912  std::unique_ptr<TFile> fTemp(TFile::Open(pythiaXSecFilename.c_str(), "READ"));
913 
914  if (!fTemp) {
915  AliDebugStream(4) << "File " << pythiaXSecFilename.c_str() << " does not exist!\n";
916  pythiaXSecFilename = "";
917  }
918  else {
919  AliDebugStream(4) << "Found pythia cross section file \"" << pythiaXSecFilename.c_str() << "\".\n";
920  }
921  }
922  }
923  else
924  {
925  // Hanlde normal root files
926  pythiaXSecFilename = gSystem->DirName(baseFileName);
927  pythiaXSecFilename += "/";
928  pythiaXSecFilename += pythiaBaseFilename;
929 
930  // Check if the file is accessible
931  if (testIfExists) {
932  if (gSystem->AccessPathName(pythiaXSecFilename.c_str())) {
933  AliDebugStream(4) << "File " << pythiaXSecFilename.c_str() << " does not exist!\n";
934  pythiaXSecFilename = "";
935  }
936  else {
937  AliDebugStream(4) << "Found pythia cross section file \"" << pythiaXSecFilename.c_str() << "\".\n";
938  }
939  }
940  }
941 
942  return pythiaXSecFilename;
943 }
944 
951 {
952  if (fInitializedConfiguration == false) {
953  AliFatal("The configuration is not initialized. Check that Initialize() was called!");
954  }
955 
956  // Setup TChain
957  Bool_t res = SetupInputFiles();
958  if (!res) { return; }
959 
960  // Note if getting random event access
962  AliInfo("Random event number access enabled!");
963  }
964 
965  fInitializedEmbedding = kTRUE;
966 }
967 
979 {
980  // Load first entry of the (next) file so that we can query information about it
981  // (it is unaccessible otherwise).
982  // Since fUpperEntry is the total number of entries, loading it will retrieve the
983  // next tree (in the next file) since entries are indexed starting from 0.
984  fChain->GetEntry(fUpperEntry);
985 
986  // Determine tree size and current entry
987  // Set the limits of the new tree
989  // Fine to be += as long as we started at 0
990  fUpperEntry += fChain->GetTree()->GetEntries();
991 
992  // Jump ahead at random if desired
993  // Determines the offset into the tree
995  TRandom3 rand(0);
996  fOffset = TMath::Nint(rand.Rndm()*(fUpperEntry-fLowerEntry))-1;
997  }
998  else {
999  fOffset = 0;
1000  }
1001 
1002  // Sets which entry to start if the try
1004 
1005  // Keep track of the number of files that we have gone through
1006  // To start from 0, we only increment if fLowerEntry > 0
1007  if (fLowerEntry > 0) {
1008  fFileNumber++;
1009  }
1010 
1011  // Add to the count the number of files which were embedded
1012  fHistManager.FillTH1("fHistNumberOfFilesEmbedded", 1);
1013  fHistManager.FillTH1("fHistAbsoluteFileNumber", (fFileNumber + fFilenameIndex) % fMaxNumberOfFiles);
1014 
1015  // Check for pythia cross section and extract if possible
1016  // fFileNumber corresponds to the next file
1017  // If there are pythia filenames, the number of match the file number of the tree.
1018  // If we previously gave up on extracting then there should be no entires
1019  if (fPythiaCrossSectionFilenames.size() > 0) {
1020  // Need to check that fFileNumber is smaller than the size of the vector because we don't check if
1023 
1024  if (!success) {
1025  AliDebugStream(3) << "Failed to retrieve cross section from xsec file. Will still attempt to get the information from the header.\n";
1026  }
1027  }
1028  else {
1029  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";
1030  }
1031  }
1032 
1033  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));
1034  // NOTE: Cannot use this print message, as it is possible that fMaxNumberOfFiles != fFilenames.size() because
1035  // invalid filenames may be included in the fFilenames count!
1036  //AliDebug(2, TString::Format("Will start embedding file %i as the %ith file beginning from entry %i.", (fFilenameIndex + fFileNumber) % fMaxNumberOfFiles, fFileNumber, fCurrentEntry));
1037 
1038  // (re)set whether we have wrapped the tree
1039  fWrappedAroundTree = false;
1040 
1041  // Note that the tree in the new file has been initialized
1042  fInitializedNewFile = kTRUE;
1043 }
1044 
1053 {
1054  std::unique_ptr<TFile> fxsec(TFile::Open(pythiaFileName.c_str()));
1055 
1056  if (fxsec)
1057  {
1058  int trials = 0;
1059  double crossSection = 0;
1060  double nEvents = 0;
1061  // Check if it's a tree
1062  TTree *xtree = dynamic_cast<TTree*>(fxsec->Get("Xsection"));
1063  if (xtree) {
1064  UInt_t ntrials = 0;
1065  Double_t xsection = 0;
1066  xtree->SetBranchAddress("xsection",&xsection);
1067  xtree->SetBranchAddress("ntrials",&ntrials);
1068  xtree->GetEntry(0);
1069  trials = ntrials;
1070  crossSection = xsection;
1071  // TODO: Test this on a file which has pyxsec.root!
1072  nEvents = 1.;
1073  AliFatal("Have no tested pyxsec.root files. Need to determine the proper way to get nevents!!");
1074  }
1075  else {
1076  // Check if it's instead the histograms
1077  // find the tlist we want to be independtent of the name so use the Tkey
1078  TKey* key = static_cast<TKey*>(fxsec->GetListOfKeys()->At(0));
1079  if (!key) return false;
1080  TList *list = dynamic_cast<TList*>(key->ReadObj());
1081  if (!list) return false;
1082  TProfile * crossSectionHist = static_cast<TProfile*>(list->FindObject("h1Xsec"));
1083  // check for failure
1084  if(!(crossSectionHist->GetEntries())) {
1085  // No cross seciton information available - fall back to raw
1086  AliErrorStream() << "No cross section information available in file \"" << fxsec->GetName() << "\". Will still attempt to extract cross section information from pythia header.\n";
1087  } else {
1088  // Cross section histogram filled - take it from there
1089  crossSection = crossSectionHist->GetBinContent(1);
1090  if(!crossSection) AliErrorStream() << GetName() << ": Cross section 0 for file " << pythiaFileName << std::endl;
1091  }
1092  TH1F * trialsHist = static_cast<TH1F*>(list->FindObject("h1Trials"));
1093  trials = trialsHist->GetBinContent(1);
1094  nEvents = trialsHist->GetEntries();
1095  }
1096 
1097  // If successful in retrieveing the values, normalizae the xsec and trials by the number of events
1098  // in the file. This way, we can use it as an approximate event-by-event value
1099  // We do not want to just use the overall value because some of the events may be rejected by various
1100  // event selections, so we only want that ones that were actually use. The easiest way to do so is by
1101  // filling it for each event.
1102  fPythiaTrialsFromFile = trials/nEvents;
1103  // Do __NOT__ divide by nEvents here! The value is already from a TProfile and therefore is already the mean!
1104  fPythiaCrossSectionFromFile = crossSection;
1105 
1106  return true;
1107  }
1108  else {
1109  AliDebugStream(3) << "Unable to open file \"" << pythiaFileName << "\". Will attempt to use values from the hader.";
1110  }
1111 
1112  // Could not open file
1113  return false;
1114 }
1115 
1122 {
1123  if (!fInitializedEmbedding) {
1124  AliError("Chain not initialized before running! Setting up now.");
1125  SetupEmbedding();
1126  }
1127 
1128  if (!fInitializedNewFile) {
1129  InitTree();
1130  }
1131 
1132  Bool_t res = GetNextEntry();
1133 
1134  if (!res) {
1135  AliError("Unable to get the event to embed. Nothing will be embedded.");
1136  return;
1137  }
1138 
1139  if (fCreateHisto && fOutput) {
1140  PostData(1, fOutput);
1141  }
1142 }
1143 
1148 {
1149 }
1150 
1159 {
1160  // Get the pointer to the existing analysis manager via the static access method.
1161  //==============================================================================
1162  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1163  if (!mgr)
1164  {
1165  ::Error("AddTaskEmcalEmbeddingHelper", "No analysis manager to connect to.");
1166  return 0;
1167  }
1168 
1169  // Check the analysis type using the event handlers connected to the analysis manager.
1170  //==============================================================================
1171  AliVEventHandler* handler = mgr->GetInputEventHandler();
1172  if (!handler)
1173  {
1174  ::Error("AddTaskEmcalEmbeddingHelper", "This task requires an input event handler");
1175  return 0;
1176  }
1177 
1178  TString name = "AliAnalysisTaskEmcalEmbeddingHelper";
1179 
1180  AliAnalysisTaskEmcalEmbeddingHelper * mgrTask = static_cast<AliAnalysisTaskEmcalEmbeddingHelper *>(mgr->GetTask(name.Data()));
1181  if (mgrTask) return mgrTask;
1182 
1183  // Create the task that manages
1184  AliAnalysisTaskEmcalEmbeddingHelper * embeddingHelper = new AliAnalysisTaskEmcalEmbeddingHelper(name.Data());
1185 
1186  //-------------------------------------------------------
1187  // Final settings, pass to manager and set the containers
1188  //-------------------------------------------------------
1189 
1190  mgr->AddTask(embeddingHelper);
1191 
1192  // Create containers for input/output
1193  AliAnalysisDataContainer* cInput = mgr->GetCommonInputContainer();
1194 
1195  TString outputContainerName(name);
1196  outputContainerName += "_histos";
1197 
1198  AliAnalysisDataContainer * cOutput = mgr->CreateContainer(outputContainerName.Data(),
1199  TList::Class(),
1200  AliAnalysisManager::kOutputContainer,
1201  Form("%s", AliAnalysisManager::GetCommonFileName()));
1202 
1203  mgr->ConnectInput(embeddingHelper, 0, cInput);
1204  mgr->ConnectOutput(embeddingHelper, 1, cOutput);
1205 
1206  return embeddingHelper;
1207 }
1208 
1214 std::string AliAnalysisTaskEmcalEmbeddingHelper::toString(bool includeFileList) const
1215 {
1216  std::stringstream tempSS;
1217 
1218  // Show the correction components
1219  tempSS << std::boolalpha;
1220  tempSS << GetName() << ": Embedding helper configuration:\n";
1221  tempSS << "Create histos: " << fCreateHisto << "\n";
1222  tempSS << "Pt Hard Bin: " << fPtHardBin << "\n";
1223  tempSS << "N Pt Hard Bins: " << fNPtHardBins << "\n";
1224  tempSS << "Anchor Run: " << fAnchorRun << "\n";
1225  tempSS << "File pattern: \"" << fFilePattern << "\"\n";
1226  tempSS << "Input filename: \"" << fInputFilename << "\"\n";
1227  tempSS << "File list filename: \"" << fFileListFilename << "\"\n";
1228  tempSS << "Tree name: " << fTreeName << "\n";
1229  tempSS << "Random event number access: " << fRandomEventNumberAccess << "\n";
1230  tempSS << "Random file access: " << fRandomFileAccess << "\n";
1231  tempSS << "Starting file index: " << fFilenameIndex << "\n";
1232  tempSS << "Number of files to embed: " << fFilenames.size() << "\n";
1233 
1234  std::bitset<32> triggerMask(fTriggerMask);
1235  tempSS << "\nEmbedded event settings:\n";
1236  tempSS << "Trigger mask (binary): " << triggerMask << "\n";
1237  tempSS << "Reject outliers: " << fMCRejectOutliers << "\n";
1238  tempSS << "Pt hard jet pt rejection factor: " << fPtHardJetPtRejectionFactor << "\n";
1239  tempSS << "Z vertex cut: " << fZVertexCut << "\n";
1240  tempSS << "Max vertex distance: " << fMaxVertexDist << "\n";
1241 
1242  if (includeFileList) {
1243  tempSS << "\nFiles to embed:\n";
1244  for (auto filename : fFilenames) {
1245  tempSS << "\t" << filename << "\n";
1246  }
1247  }
1248 
1249  return tempSS.str();
1250 }
1251 
1259 std::ostream & AliAnalysisTaskEmcalEmbeddingHelper::Print(std::ostream & in) const {
1260  in << toString();
1261  return in;
1262 }
1263 
1272 std::ostream & operator<<(std::ostream & in, const AliAnalysisTaskEmcalEmbeddingHelper & myTask)
1273 {
1274  std::ostream & result = myTask.Print(in);
1275  return result;
1276 }
1277 
1285 {
1286  std::string temp(opt);
1287  bool includeFileList = false;
1288  if (temp == "FILELIST") {
1289  includeFileList = true;
1290  }
1291  Printf("%s", toString(includeFileList).c_str());
1292 }
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