AliPhysics  cdeda5a (cdeda5a)
 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 <vector>
18 #include <algorithm>
19 #include <fstream>
20 
21 #include <TFile.h>
22 #include <TMath.h>
23 #include <TRandom.h>
24 #include <TChain.h>
25 #include <TGrid.h>
26 #include <TSystem.h>
27 #include <TUUID.h>
28 
29 #include <AliLog.h>
30 #include <AliAnalysisManager.h>
31 #include <AliVEvent.h>
32 #include <AliAODEvent.h>
33 #include <AliESDEvent.h>
34 #include <AliInputEventHandler.h>
35 
37 
41 
43 
49  fTreeName(),
50  fAnchorRun(169838),
51  fPtHardBin(-1),
52  fRandomEventNumberAccess(kFALSE),
53  fRandomFileAccess(kFALSE),
54  fFilePattern(""),
55  fFileListFilename(""),
56  fFilenameIndex(-1),
57  fFilenames(),
58  fTriggerMask(AliVEvent::kAny),
59  fZVertexCut(10),
60  fMaxVertexDist(999),
61  fExternalFile(0),
62  fCurrentEntry(0),
63  fLowerEntry(0),
64  fUpperEntry(0),
65  fOffset(0),
66  fMaxNumberOfFiles(0),
67  fFileNumber(0),
68  fInitializedEmbedding(false),
69  fInitializedNewFile(false),
70  fWrappedAroundTree(false),
71  fChain(0),
72  fExternalEvent(0)
73 {
74  if (fgInstance != 0) {
75  AliError("An instance of AliAnalysisTaskEmcalEmbeddingHelper already exists: it will be deleted!!!");
76  delete fgInstance;
77  }
78 
79  fgInstance = this;
80 }
81 
88  AliAnalysisTaskSE(name),
89  fTreeName("aodTree"),
90  fAnchorRun(169838),
91  fPtHardBin(-1),
92  fRandomEventNumberAccess(kFALSE),
93  fRandomFileAccess(kFALSE),
94  fFilePattern(""),
95  fFileListFilename(""),
96  fFilenameIndex(-1),
97  fFilenames(),
98  fTriggerMask(AliVEvent::kAny),
99  fZVertexCut(10),
100  fMaxVertexDist(999),
101  fExternalFile(0),
102  fCurrentEntry(0),
103  fLowerEntry(0),
104  fUpperEntry(0),
105  fOffset(0),
106  fMaxNumberOfFiles(0),
107  fFileNumber(0),
108  fInitializedEmbedding(false),
109  fInitializedNewFile(false),
110  fWrappedAroundTree(false),
111  fChain(0),
112  fExternalEvent(0)
113 {
114  if (fgInstance != 0) {
115  AliError("An instance of AliAnalysisTaskEmcalEmbeddingHelper already exists: it will be deleted!!!");
116  delete fgInstance;
117  }
118 
119  fgInstance = this;
120 }
121 
128 {
129  if (fgInstance == this) fgInstance = 0;
130  if (fExternalEvent) delete fExternalEvent;
131  if (fExternalFile) {
132  fExternalFile->Close();
133  delete fExternalFile;
134  }
135 }
136 
161 {
162  // Retrieve filenames if we don't have them yet.
163  if (fFilenames.size() == 0)
164  {
165  // Handle if fPtHardBin or fAnchorRun are set
166  // This will require formatting the file pattern in the proper way to support these substitutions
167  if (fPtHardBin != -1) {
168  if (fAnchorRun > 0) {
169  fFilePattern = TString::Format(fFilePattern, fAnchorRun, fPtHardBin);
170  }
171  else {
172  fFilePattern = TString::Format(fFilePattern, fPtHardBin);
173  }
174  }
175 
176  // Retrieve AliEn filenames
177  if (fFilePattern.BeginsWith("alien://") && !gGrid) {
178  AliInfo("Trying to connect to AliEn ...");
179  TGrid::Connect("alien://");
180 
181  if (!gGrid) {
182  AliFatal(Form("Cannot access AliEn to retrieve file list with pattern %s!", fFilePattern.Data()));
183  }
184 
185  AliDebug(2,Form("Trying to retrieve file list from AliEn with pattern file %s...", fFilePattern.Data()));
186 
187  // Create a temporary filename based on a UUID to make sure that it doesn't overwrite any files
188  if (fFileListFilename == "") {
189  TUUID tempUUID;
190  fFileListFilename = "fileList";
191  fFileListFilename += tempUUID.AsString();
192  fFileListFilename += ".txt";
193  }
194 
195  // Retrieve filenames from alien using alien_find
196  TString command = "alien_find -l";
197  command += fFilePattern;
198  command += " > ";
199  command += fFileListFilename;
200 
201  // Execute the alien_find command to get the filenames
202  AliDebug(2,Form("Trying to retrieve file list from AliEn with alien_find command \"%s\"", command.Data()));
203  gSystem->Exec(command.Data());
204  }
205 
206  std::ifstream inputFile(fFileListFilename);
207 
208  // Copy available files into the filenames vector
209  // From:: https://stackoverflow.com/a/8365247
210  std::copy(std::istream_iterator<std::string>(inputFile),
211  std::istream_iterator<std::string>(),
212  std::back_inserter(fFilenames));
213  }
214 
215  if (fFilenames.size() == 0) {
216  AliFatal(Form("Filenames from pattern \"%s\" and file list \"%s\" yielded an empty list!", fFilePattern.Data(), fFileListFilename.Data()));
217  }
218 
219  // Add "#" to files in there are any zip files
220  // It is require to open the proper root file within the zip
221  for (auto filename : fFilenames)
222  {
223  if (filename.find(".zip") != std::string::npos && filename.find("#") == std::string::npos) {
224  if (fTreeName == "aodTree") {
225  filename += "#AliAOD.root";
226  }
227  else if (fTreeName == "esdTree") {
228  filename += "#AliESDs.root";
229  }
230  else {
231  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()));
232  }
233  }
234  }
235 
236  // Determine the next filename index
237  // This determines which file is added first to the TChain, thus determining the order of processing
238  // Random file access. Only do this if the user has no set the filename index and request random file access
239  if (fFilenameIndex == -1 && fRandomFileAccess) {
240  // - 1 ensures that we it doesn't overflow
241  fFilenameIndex = TMath::Nint(gRandom->Rndm()*fFilenames.size()) - 1;
242  AliInfo(TString::Format("Starting with random file number %i!", fFilenameIndex));
243  }
244  // If not random file access, then start from the beginning
245  if (fFilenameIndex >= fFilenames.size() || fFilenameIndex < 0) {
246  // Skip notifying on -1 since it will likely be set there due to constructor.
247  if (fFilenameIndex != -1) {
248  AliWarning(Form("File index %i out of range from 0 to %lu! Resetting to 0!", fFilenameIndex, fFilenames.size()));
249  }
250  fFilenameIndex = 0;
251  }
252 
253  AliInfo(TString::Format("Starting with file number %i out of %lu", fFilenameIndex, fFilenames.size()));
254 }
255 
265 {
266  Int_t attempts = -1;
267 
268  do {
269  // Reset to start of tree
270  if (fCurrentEntry == fUpperEntry) {
272  fWrappedAroundTree = true;
273  }
274 
276  // Continue with GetEntry as normal
277  }
278  else {
279  // NOTE: On transition from one file to the next, this calls the next entry that would be expected.
280  // However, if it is for the last file, it tries to GetEntry() of one entry past the end of the last file.
281  // Normally, this would be a problem, however GetEntry() just doesn't fill the fields of an invalid index
282  // instead of throwing an error. So "invalid values" are filled for a file that doesn't exist, but then
283  // they are immediately replaced by the lines below that reset the access values and re-init the tree.
284  // The benefit of this approach is it simplies file counting (we don't need to carefully increment here
285  // and in InitTree()) and preserves the desired behavior when we are not at the last file.
286  InitTree();
287  }
288 
289  // Load current event
290  // Can be a simple less than, because fFileNumber counts from 0.
292  fChain->GetEntry(fCurrentEntry);
293  }
294  else {
295  AliError("====================================================================================================");
296  AliError("== No more files available to embed from the TChain! Restarting from the beginning of the TChain! ==");
297  AliError("== Be careful to check that this is the desired action! ==");
298  AliError("====================================================================================================");
299 
300  // Reset the relevant access values
301  // fCurrentEntry and fLowerEntry are automatically reset in InitTree()
302  fFileNumber = 0;
303  fUpperEntry = 0;
304 
305  // Re-init back to the start
306  InitTree();
307 
308  // Access the relevant entry
309  // We are certain that fFileNumber is less than fMaxNumberOfFiles, so we are resetting to start
310  fChain->GetEntry(fCurrentEntry);
311  }
312  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));
313 
314  // Increment current entry
315  fCurrentEntry++;
316 
317  // Provide a check for number of attempts
318  attempts++;
319  if (attempts == 1000)
320  AliWarning("After 1000 attempts no event has been accepted by the event selection (trigger, centrality...)!");
321 
322  } while (!IsEventSelected());
323 
324  if (!fChain) return kFALSE;
325 
326  return kTRUE;
327 }
328 
335 {
336  // AOD Event selection.
337  // TODO: Can we make centrality make sense here?
338  /*if (!fEsdTreeMode && fAODHeader) {
339  AliAODHeader *aodHeader = static_cast<AliAODHeader*>(fAODHeader);
340 
341  // Centrality selection
342  if (fMinCentrality >= 0) {
343  AliCentrality *cent = aodHeader->GetCentralityP();
344  Float_t centVal = cent->GetCentralityPercentile("V0M");
345  if (centVal < fMinCentrality || centVal >= fMaxCentrality) {
346  AliDebug(2,Form("Event rejected due to centrality selection. Event centrality: %f, centrality range selection: %f to %f",
347  centVal, fMinCentrality, fMaxCentrality));
348  return kFALSE;
349  }
350  }
351  }*/
352 
353  // Trigger selection
354  if (fTriggerMask != AliVEvent::kAny) {
355  UInt_t res = 0;
356  const AliESDEvent *eev = dynamic_cast<const AliESDEvent*>(InputEvent());
357  if (eev) {
358  res = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
359  } else {
360  const AliAODEvent *aev = dynamic_cast<const AliAODEvent*>(InputEvent());
361  if (aev) {
362  res = ((AliVAODHeader*)aev->GetHeader())->GetOfflineTrigger();
363  }
364  }
365 
366  if ((res & fTriggerMask) == 0) {
367  AliDebug(3, Form("Event rejected due to physics selection. Event trigger mask: %d, trigger mask selection: %d.",
368  res, fTriggerMask));
369  return kFALSE;
370  }
371  }
372 
373  // Vertex selection
374  Double_t externalVertex[3]={0};
375  Double_t inputVertex[3]={0};
376  const AliVVertex *externalVert = fExternalEvent->GetPrimaryVertex();
377  const AliVVertex *inputVert = InputEvent()->GetPrimaryVertex();
378  if (externalVert && inputVert) {
379  externalVert->GetXYZ(externalVertex);
380  inputVert->GetXYZ(inputVertex);
381 
382  if (TMath::Abs(externalVertex[2]) > fZVertexCut) {
383  AliDebug(3, Form("Event rejected due to Z vertex selection. Event Z vertex: %f, Z vertex cut: %f",
384  externalVertex[2], fZVertexCut));
385  return kFALSE;
386  }
387  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]));
388  if (dist > fMaxVertexDist) {
389  AliDebug(3, Form("Event rejected because the distance between the current and embedded vertices is > %f. "
390  "Current event vertex (%f, %f, %f), embedded event vertex (%f, %f, %f). Distance = %f",
391  fMaxVertexDist, inputVertex[0], inputVertex[1], inputVertex[2], externalVertex[0], externalVertex[1], externalVertex[2], dist));
392  return kFALSE;
393  }
394  }
395 
396  // TODO: Can we do selection based on the contents of the external event input objects?
397  // The previous embedding task could do so by directly accessing the elements.
398  // Certainly can't do jets (say minPt of leading jet) :because this has to be embedded before them.
399  // See AliJetEmbeddingFromAODTask::IsAODEventSelected()
400 
401  return kTRUE;
402 }
403 
410 {
411  if (!fChain) return kFALSE;
412 
413  if (!fExternalEvent) {
414  if (fTreeName == "aodTree") {
415  fExternalEvent = new AliAODEvent();
416  }
417  else if (fTreeName == "esdTree") {
418  fExternalEvent = new AliESDEvent();
419  }
420  else {
421  AliError(Form("Tree name %s not recognized!", fTreeName.Data()));
422  return kFALSE;
423  }
424  }
425 
426  fExternalEvent->ReadFromTree(fChain, fTreeName);
427 
428  return kTRUE;
429 }
430 
437 {
438  SetupEmbedding();
439 }
440 
448 {
449  // Setup TChain
450  fChain = new TChain(fTreeName);
451 
452  // Determine whether AliEn is needed
453  bool requiresAlien = false;
454  for (auto filename : fFilenames)
455  {
456  if (filename.find("alien://") != std::string::npos) {
457  requiresAlien = true;
458  }
459  }
460 
461  if (requiresAlien && !gGrid) {
462  AliInfo("Trying to connect to AliEn ...");
463  TGrid::Connect("alien://");
464  }
465 
466  // Add files for TChain
467  // See: https://stackoverflow.com/a/8533198
468  bool wrapped = false;
469  TString baseFileName("");
470  for (auto filename = fFilenames.begin() + fFilenameIndex; (filename != fFilenames.begin() + fFilenameIndex || !wrapped); filename++)
471  {
472  // Wraps the loop back around to the beginning
473  if (filename == fFilenames.end())
474  {
475  // Explicit check is needed. Otherwise, an offset of 0 would load the 0th entry twice.
476  if (fFilenameIndex == 0) {
477  break;
478  }
479  filename = fFilenames.begin();
480  wrapped = true;
481  }
482 
483  // AccessPathName() cannot handle the "#", so we need to strip it to check that the file exists.
484  baseFileName = filename->c_str();
485  if (baseFileName.Contains(".zip#")) {
486  Ssiz_t pos = baseFileName.Last('#');
487  baseFileName.Remove(pos);
488  }
489 
490  // Ensure that the file is accessible
491  if (gSystem->AccessPathName(baseFileName)) {
492  AliError(Form("File %s does not exist! Skipping!", baseFileName.Data()));
493  // Do not process the file if it is unaccessible, but continue processing
494  continue;
495  }
496 
497  // Add to the Chain
498  fChain->Add(filename->c_str());
499  }
500 
501  // Keep track of the total number of files in the TChain to ensure that we don't start repeating within the chain
502  fMaxNumberOfFiles = fChain->GetListOfFiles()->GetEntries();
503 
504  if (fFilenames.size() > fMaxNumberOfFiles) {
505  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));
506  }
507 
508  // Setup input event
509  Bool_t res = InitEvent();
510  if (!res) return kFALSE;
511 
512  return kTRUE;
513 }
514 
521 {
522  // Get file list
523  GetFilenames();
524 
525  // Setup TChain
526  Bool_t res = SetupInputFiles();
527  if (!res) { return; }
528 
529  // Note if getting random event access
531  AliInfo("Random event number access enabled!");
532  }
533 
534  fInitializedEmbedding = kTRUE;
535 }
536 
548 {
549  // Load first entry of the (next) file so that we can query information about it
550  // (it is unaccessible otherwise).
551  // Since fUpperEntry is the total number of entries, loading it will retrieve the
552  // next tree (in the next file) since entries are indexed starting from 0.
553  fChain->GetEntry(fUpperEntry);
554 
555  // Determine tree size and current entry
556  // Set the limits of the new tree
558  // Fine to be += as long as we started at 0
559  fUpperEntry += fChain->GetTree()->GetEntries();
560 
561  // Jump ahead at random if desired
562  // Determines the offset into the tree
564  fOffset = TMath::Nint(gRandom->Rndm()*(fUpperEntry-fLowerEntry))-1;
565  }
566  else {
567  fOffset = 0;
568  }
569 
570  // Sets which entry to start if the try
572 
573  // Keep track of the number of files that we have gone through
574  // To start from 0, we only increment if fLowerEntry > 0
575  if (fLowerEntry > 0) {
576  fFileNumber++;
577  }
578 
579  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));
580  // NOTE: Cannot use this print message, as it is possible that fMaxNumberOfFiles != fFilenames.size() because
581  // invalid filenames may be included in the fFilenames count!
582  //AliDebug(2, TString::Format("Will start embedding file %i as the %ith file beginning from entry %i.", (fFilenameIndex + fFileNumber) % fMaxNumberOfFiles, fFileNumber, fCurrentEntry));
583 
584  // (re)set whether we have wrapped the tree
585  fWrappedAroundTree = false;
586 
587  // Note that the tree in the new file has been initialized
588  fInitializedNewFile = kTRUE;
589 }
590 
597 {
598  if (!fInitializedEmbedding) {
599  AliError("Chain not initialized before running! Setting up now.");
600  SetupEmbedding();
601  }
602 
603  if (!fInitializedNewFile) {
604  InitTree();
605  }
606 
607  Bool_t res = GetNextEntry();
608 
609  if (!res) {
610  AliError("Unable to get the event to embed. Nothing will be embedded.");
611  return;
612  }
613 }
614 
619 {
620 }
621 
630 {
631  // Get the pointer to the existing analysis manager via the static access method.
632  //==============================================================================
633  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
634  if (!mgr)
635  {
636  ::Error("AddTaskEmcalEmbeddingHelper", "No analysis manager to connect to.");
637  return 0;
638  }
639 
640  // Check the analysis type using the event handlers connected to the analysis manager.
641  //==============================================================================
642  AliVEventHandler* handler = mgr->GetInputEventHandler();
643  if (!handler)
644  {
645  ::Error("AddTaskEmcalEmbeddingHelper", "This task requires an input event handler");
646  return 0;
647  }
648 
649  TString name = "AliAnalysisTaskEmcalEmbeddingHelper";
650 
651  AliAnalysisTaskEmcalEmbeddingHelper * mgrTask = static_cast<AliAnalysisTaskEmcalEmbeddingHelper *>(mgr->GetTask(name.Data()));
652  if (mgrTask) return mgrTask;
653 
654  // Create the task that manages
656 
657  //-------------------------------------------------------
658  // Final settings, pass to manager and set the containers
659  //-------------------------------------------------------
660 
661  mgr->AddTask(embeddingHelper);
662 
663  // Create containers for input/output
664  AliAnalysisDataContainer* cInput = mgr->GetCommonInputContainer();
665 
666  /*TString outputContainerName(name);
667  outputContainerName += "_histos";
668 
669  AliAnalysisDataContainer * cOutput = mgr->CreateContainer(outputContainerName.Data(),
670  TList::Class(),
671  AliAnalysisManager::kOutputContainer,
672  Form("%s", AliAnalysisManager::GetCommonFileName()));*/
673 
674  mgr->ConnectInput(embeddingHelper, 0, cInput);
675  //mgr->ConnectOutput(embeddingHelper, 1, cOutput);
676 
677  return embeddingHelper;
678 }
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
Int_t fLowerEntry
! First entry of the current tree to be used for embedding
Int_t fAnchorRun
Anchor run for the given pythia production.
TString fFileListFilename
Name of the file list containing paths to files to embed.
Bool_t fRandomFileAccess
If true, it will start embedding from a random file in the input files list.
TSystem * gSystem
Int_t fCurrentEntry
! Current entry in the current tree
TFile * fExternalFile
! External file used for embedding
TChain * fChain
! External TChain (tree) containing the events available for embedding
Declaration of class AliAnalysisTaskEmcalEmbeddingHelper.
TRandom * gRandom
Int_t fMaxNumberOfFiles
! Max number of files that are in the TChain
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
bool fInitializedEmbedding
! Notes where the TChain has been initialized for embedding
Implementation of task to embed external events.
Double_t fMaxVertexDist
Max distance between Z vertex of internal and embedded event.
AliVEvent * fExternalEvent
! Current external event available for embedding
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 fWrappedAroundTree
! Notes whether we have wrapped around the tree, which is important if the offset into the tree is no...
static AliAnalysisTaskEmcalEmbeddingHelper * AddTaskEmcalEmbeddingHelper()
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
bool fInitializedNewFile
! Notes where the entry indices have been initialized for a new tree in the chain ...
static AliAnalysisTaskEmcalEmbeddingHelper * fgInstance
! Global instance of this class
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53
Int_t fUpperEntry
! Last entry of the current tree to be used for embedding
Int_t fFilenameIndex
Index of vector containing paths to files to embed.
std::vector< std::string > fFilenames
! Paths to the files to embed
TString fFilePattern
File pattern to select AliEn files.
Int_t fFileNumber
! File number corresponding to the current tree