AliPhysics  master (3d17d9d)
AliAnalysisTaskEmcalLight.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 #include <sstream>
16 #include <array>
17 #include <memory>
18 
19 #include <RVersion.h>
20 #include <TClonesArray.h>
21 #include <TList.h>
22 #include <TObject.h>
23 #include <TH1F.h>
24 #include <TH2F.h>
25 #include <TProfile.h>
26 #include <TSystem.h>
27 #include <TFile.h>
28 #include <TChain.h>
29 #include <TKey.h>
30 
31 #include "AliGenCocktailEventHeader.h"
32 #include "AliStack.h"
33 #include "AliAODEvent.h"
34 #include "AliAnalysisManager.h"
35 #include "AliCentrality.h"
36 #include "AliEmcalList.h"
37 #include "AliEMCALGeometry.h"
38 #include "AliESDEvent.h"
39 #include "AliEmcalParticle.h"
40 #include "AliEventplane.h"
41 #include "AliInputEventHandler.h"
42 #include "AliLog.h"
43 #include "AliMCParticle.h"
44 #include "AliVCluster.h"
45 #include "AliVEventHandler.h"
46 #include "AliVParticle.h"
47 #include "AliAODTrack.h"
48 #include "AliVCaloTrigger.h"
49 #include "AliGenPythiaEventHeader.h"
50 #include "AliGenEventHeader.h"
51 #include "AliAODMCHeader.h"
52 #include "AliMCEvent.h"
53 #include "AliEMCALTriggerPatchInfo.h"
54 
55 #include "AliMultSelection.h"
56 
58 
60 
62 
65  fForceBeamType(kNA),
66  fGeneralHistograms(kFALSE),
67  fCreateHisto(kTRUE),
68  fNeedEmcalGeom(kTRUE),
69  fUseBuiltinEventSelection(kFALSE),
70  fCentBins(),
71  fCentralityEstimation(kNewCentrality),
72  fIsPythia(kFALSE),
73  fIsMonteCarlo(kFALSE),
74  fMCEventHeaderName(),
75  fCaloCellsName(),
76  fCaloTriggersName(),
77  fCaloTriggerPatchInfoName(),
78  fCentEst("V0M"),
79  fParticleCollArray(),
80  fClusterCollArray(),
81  fTriggerSelectionBitMap(0),
82  fMinCent(-1),
83  fMaxCent(-1),
84  fMinVz(-999),
85  fMaxVz(999),
86  fMaxVzDiff(-1),
87  fMinNVertCont(0),
88  fMinPtHard(-1),
89  fMaxPtHard(-1),
90  fMaxMinimumBiasPtHard(-1),
91  fAcceptedTriggerClasses(),
92  fRejectedTriggerClasses(),
93  fMCRejectFilter(kFALSE),
94  fPtHardAndJetPtFactor(0.),
95  fPtHardAndClusterPtFactor(0.),
96  fPtHardAndTrackPtFactor(0.),
97  fSwitchOffLHC15oFaultyBranches(kFALSE),
98  fEventSelectionAfterRun(kFALSE),
99  fUseAliEmcalList(kFALSE),
100  fUsePtHardBinScaling(kFALSE),
101  fSelectGeneratorName(),
102  fMinimumEventWeight(1e-6),
103  fMaximumEventWeight(1e6),
104  fInhibit(kFALSE),
105  fLocalInitialized(kFALSE),
106  fWarnMissingCentrality(kTRUE),
107  fDataType(kAOD),
108  fGeom(0),
109  fCaloCells(0),
110  fCaloTriggers(0),
111  fTriggerPatchInfo(0),
112  fCent(-1),
113  fCentBin(-1),
114  fEPV0(-1.0),
115  fEPV0A(-1.0),
116  fEPV0C(-1.0),
117  fNVertCont(0),
118  fNVertSPDCont(0),
119  fFiredTriggerBitMap(0),
120  fFiredTriggerClasses(),
121  fBeamType(kNA),
122  fMCHeader(0),
123  fPythiaHeader(0),
124  fUseXsecFromHeader(false),
125  fPtHardBin(0),
126  fPtHard(0),
127  fNTrials(0),
128  fXsection(0),
129  fEventWeight(1),
130  fGeneratorName(),
131  fOutput(0),
132  fHistograms()
133 {
134  fVertex[0] = 0;
135  fVertex[1] = 0;
136  fVertex[2] = 0;
137  fVertexSPD[0] = 0;
138  fVertexSPD[1] = 0;
139  fVertexSPD[2] = 0;
140 }
141 
143  AliAnalysisTaskSE(name),
144  fForceBeamType(kNA),
145  fGeneralHistograms(kFALSE),
146  fCreateHisto(kTRUE),
147  fNeedEmcalGeom(kTRUE),
148  fUseBuiltinEventSelection(kFALSE),
149  fCentBins(6),
150  fCentralityEstimation(kNewCentrality),
151  fIsPythia(kFALSE),
152  fIsMonteCarlo(kFALSE),
153  fMCEventHeaderName(),
154  fCaloCellsName(),
155  fCaloTriggersName(),
156  fCaloTriggerPatchInfoName(),
157  fCentEst("V0M"),
158  fParticleCollArray(),
159  fClusterCollArray(),
160  fTriggerSelectionBitMap(0),
161  fMinCent(-1),
162  fMaxCent(-1),
163  fMinVz(-999),
164  fMaxVz(999),
165  fMaxVzDiff(-1),
166  fMinNVertCont(0),
167  fMinPtHard(-1),
168  fMaxPtHard(-1),
169  fMaxMinimumBiasPtHard(-1),
170  fAcceptedTriggerClasses(),
171  fRejectedTriggerClasses(),
172  fMCRejectFilter(kFALSE),
173  fPtHardAndJetPtFactor(0.),
174  fPtHardAndClusterPtFactor(0.),
175  fPtHardAndTrackPtFactor(0.),
176  fSwitchOffLHC15oFaultyBranches(kFALSE),
177  fEventSelectionAfterRun(kFALSE),
178  fUseAliEmcalList(kFALSE),
179  fUsePtHardBinScaling(kFALSE),
180  fSelectGeneratorName(),
181  fMinimumEventWeight(1e-6),
182  fMaximumEventWeight(1e6),
183  fInhibit(kFALSE),
184  fLocalInitialized(kFALSE),
185  fWarnMissingCentrality(kTRUE),
186  fDataType(kAOD),
187  fGeom(0),
188  fCaloCells(0),
189  fCaloTriggers(0),
190  fTriggerPatchInfo(0),
191  fCent(0),
192  fCentBin(-1),
193  fEPV0(-1.0),
194  fEPV0A(-1.0),
195  fEPV0C(-1.0),
196  fNVertCont(0),
197  fNVertSPDCont(0),
198  fFiredTriggerBitMap(0),
199  fFiredTriggerClasses(),
200  fBeamType(kNA),
201  fMCHeader(0),
202  fPythiaHeader(0),
203  fUseXsecFromHeader(false),
204  fPtHardBin(0),
205  fPtHard(0),
206  fNTrials(0),
207  fXsection(0),
208  fEventWeight(1),
209  fGeneratorName(),
210  fOutput(0),
211  fHistograms()
212 {
213  fVertex[0] = 0;
214  fVertex[1] = 0;
215  fVertex[2] = 0;
216  fVertexSPD[0] = 0;
217  fVertexSPD[1] = 0;
218  fVertexSPD[2] = 0;
219 
220  fCentBins[0] = 0;
221  fCentBins[1] = 10;
222  fCentBins[2] = 30;
223  fCentBins[3] = 50;
224  fCentBins[4] = 90;
225  fCentBins[5] = 100;
226 
227  fAliEventCuts.OverrideAutomaticTriggerSelection(AliVEvent::kAny, true);
228 
229  if (fCreateHisto) DefineOutput(1, TList::Class());
230 }
231 
233 {
234  for (auto cont_it : fParticleCollArray) delete cont_it.second;
235  for (auto cont_it : fClusterCollArray) delete cont_it.second;
236 }
237 
239 {
240  if (fInhibit) {
241  AliWarningStream() << "The execution of this task is inhibited. Returning." << std::endl;
242  return;
243  }
244  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
245  if (mgr) {
246  AliVEventHandler *evhand = mgr->GetInputEventHandler();
247  if (evhand) {
248  if (evhand->InheritsFrom("AliESDInputHandler") || evhand->InheritsFrom("AliDummyHandler")) {
249  fDataType = kESD;
250  }
251  else {
252  fDataType = kAOD;
253  }
254  }
255  else {
256  AliError("Event handler not found!");
257  }
258  }
259  else {
260  AliError("Analysis manager not found!");
261  }
262 
263  if (!fCreateHisto)
264  return;
265 
266  OpenFile(1);
267  if(fUseAliEmcalList) {
268  auto emclist = new AliEmcalList;
269  if(fUsePtHardBinScaling) emclist->SetUseScaling(true);
270  emclist->SetNameXsec("fHistXsectionExternalFile");
271  emclist->SetNameTrials("fHistTrialsExternalFile");
272  fOutput = emclist;
273  } else {
274  fOutput = new TList();
275  }
276  fOutput->SetOwner(); // @suppress("Ambiguous problem")
277 
279 
280  if (!fGeneralHistograms) return;
281 
282  TH1* h = nullptr;
283 
284  if (fIsMonteCarlo) {
285  auto weight_bins = GenerateLogFixedBinArray(1000, fMinimumEventWeight, fMaximumEventWeight, true);
286 
287  h = new TH1F("fHistEventsVsPtHard", "fHistEventsVsPtHard", 1000, 0, 1000);
288  h->GetXaxis()->SetTitle("#it{p}_{T,hard} (GeV/#it{c})");
289  h->GetYaxis()->SetTitle("events");
290  fOutput->Add(h);
291  fHistograms["fHistEventsVsPtHard"] = h;
292 
293  h = new TH1F("fHistTrialsVsPtHard", "fHistTrialsVsPtHard", 1000, 0, 1000);
294  h->GetXaxis()->SetTitle("#it{p}_{T,hard} (GeV/#it{c})");
295  h->GetYaxis()->SetTitle("trials");
296  fOutput->Add(h);
297  fHistograms["fHistTrialsVsPtHard"] = h;
298 
299  h = new TProfile("fHistXsection", "fHistXsection", 50, 0, 50);
300  h->GetXaxis()->SetTitle("#it{p}_{T,hard} bin");
301  h->GetYaxis()->SetTitle("total integrated cross section (mb)");
302  fOutput->Add(h);
303  fHistograms["fHistXsection"] = h;
304 
305  h = new TH1F("fHistXsectionDistribution", "fHistXsectionDistribution", 1000, &weight_bins[0]);
306  h->GetXaxis()->SetTitle("total integrated cross section (mb)");
307  h->GetYaxis()->SetTitle("events");
308  fOutput->Add(h);
309  fHistograms["fHistXsectionDistribution"] = h;
310 
311  h = new TH1F("fHistEventWeights", "fHistEventWeights", 1000, &weight_bins[0]);
312  h->GetXaxis()->SetTitle("weight");
313  h->GetYaxis()->SetTitle("events");
314  fOutput->Add(h);
315  fHistograms["fHistEventWeights"] = h;
316 
317  h = new TH2F("fHistEventWeightsVsPtHard", "fHistEventWeightsVsPtHard", 1000, 0, 1000, 1000, &weight_bins[0]);
318  h->GetXaxis()->SetTitle("#it{p}_{T,hard} (GeV/#it{c})");
319  h->GetYaxis()->SetTitle("event weight");
320  fOutput->Add(h);
321  fHistograms["fHistEventWeightsVsPtHard"] = h;
322 
323  h = new TH1F("fHistEventsVsPtHardNoSel", "fHistEventsVsPtHardNoSel", 1000, 0, 1000);
324  h->GetXaxis()->SetTitle("#it{p}_{T,hard} (GeV/#it{c})");
325  h->GetYaxis()->SetTitle("events");
326  fOutput->Add(h);
327  fHistograms["fHistEventsVsPtHardNoSel"] = h;
328 
329  h = new TH1F("fHistTrialsVsPtHardNoSel", "fHistTrialsVsPtHardNoSel", 1000, 0, 1000);
330  h->GetXaxis()->SetTitle("#it{p}_{T,hard} (GeV/#it{c})");
331  h->GetYaxis()->SetTitle("trials");
332  fOutput->Add(h);
333  fHistograms["fHistTrialsVsPtHardNoSel"] = h;
334 
335  h = new TProfile("fHistXsectionNoSel", "fHistXsectionNoSel", 50, 0, 50);
336  h->GetXaxis()->SetTitle("#it{p}_{T,hard} bin");
337  h->GetYaxis()->SetTitle("total integrated cross section (mb)");
338  fOutput->Add(h);
339  fHistograms["fHistXsectionNoSel"] = h;
340 
341  h = new TH1F("fHistXsectionDistributionNoSel", "fHistXsectionDistributionNoSel", 1000, &weight_bins[0]);
342  h->GetXaxis()->SetTitle("total integrated cross section (mb)");
343  h->GetYaxis()->SetTitle("events");
344  fOutput->Add(h);
345  fHistograms["fHistXsectionDistributionNoSel"] = h;
346 
347  h = new TH1F("fHistEventWeightsNoSel", "fHistEventWeightsNoSel", 1000, &weight_bins[0]);
348  h->GetXaxis()->SetTitle("weight");
349  h->GetYaxis()->SetTitle("events");
350  fOutput->Add(h);
351  fHistograms["fHistEventWeightsNoSel"] = h;
352 
353  h = new TH2F("fHistEventWeightsVsPtHardNoSel", "fHistEventWeightsVsPtHardNoSel", 1000, 0, 1000, 1000, &weight_bins[0]);
354  h->GetXaxis()->SetTitle("#it{p}_{T,hard} (GeV/#it{c})");
355  h->GetYaxis()->SetTitle("event weight");
356  fOutput->Add(h);
357  fHistograms["fHistEventWeightsVsPtHardNoSel"] = h;
358 
359  h = new TH1F("fHistTrialsExternalFile", "fHistTrialsExternalFile", 50, 0, 50);
360  h->GetXaxis()->SetTitle("#it{p}_{T,hard} bin");
361  h->GetYaxis()->SetTitle("trials");
362  fOutput->Add(h);
363  fHistograms["fHistTrialsExternalFile"] = h;
364 
365  h = new TH1F("fHistEventsExternalFile", "fHistEventsExternalFile", 50, 0, 50);
366  h->GetXaxis()->SetTitle("#it{p}_{T,hard} bin");
367  h->GetYaxis()->SetTitle("total events");
368  fOutput->Add(h);
369  fHistograms["fHistEventsExternalFile"] = h;
370 
371  h = new TProfile("fHistXsectionExternalFile", "fHistXsectionExternalFile", 50, 0, 50);
372  h->GetXaxis()->SetTitle("#it{p}_{T,hard} bin");
373  h->GetYaxis()->SetTitle("total integrated cross section (mb)");
374  fOutput->Add(h);
375  fHistograms["fHistXsectionExternalFile"] = h;
376  }
377 
378  h = new TH1F("fHistZVertex","Z vertex position", 60, -30, 30);
379  h->GetXaxis()->SetTitle("V_{#it{z}}");
380  h->GetYaxis()->SetTitle("counts");
381  fOutput->Add(h);
382  fHistograms["fHistZVertex"] = h;
383 
384  h = new TH1F("fHistZVertexNoSel","Z vertex position (no event selection)", 60, -30, 30);
385  h->GetXaxis()->SetTitle("V_{#it{z}}");
386  h->GetYaxis()->SetTitle("counts");
387  fOutput->Add(h);
388  fHistograms["fHistZVertexNoSel"] = h;
389 
391  h = new TH1F("fHistCentrality","Event centrality distribution", 100, 0, 100);
392  h->GetXaxis()->SetTitle("Centrality (%)");
393  h->GetYaxis()->SetTitle("counts");
394  fOutput->Add(h);
395  fHistograms["fHistCentrality"] = h;
396 
397  h = new TH1F("fHistCentralityNoSel","Event centrality distribution (no event selection)", 100, 0, 100);
398  h->GetXaxis()->SetTitle("Centrality (%)");
399  h->GetYaxis()->SetTitle("counts");
400  fOutput->Add(h);
401  fHistograms["fHistCentralityNoSel"] = h;
402  }
403 
404  if (fForceBeamType != kpp) {
405  h = new TH1F("fHistEventPlane","Event plane", 120, -TMath::Pi(), TMath::Pi());
406  h->GetXaxis()->SetTitle("event plane");
407  h->GetYaxis()->SetTitle("counts");
408  fOutput->Add(h);
409  fHistograms["fHistEventPlane"] = h;
410 
411  h = new TH1F("fHistEventPlaneNoSel","Event plane (no event selection)", 120, -TMath::Pi(), TMath::Pi());
412  h->GetXaxis()->SetTitle("event plane");
413  h->GetYaxis()->SetTitle("counts");
414  fOutput->Add(h);
415  fHistograms["fHistEventPlaneNoSel"] = h;
416  }
417 
418  h = new TH1F("fHistEventRejection","Reasons to reject event",30,0,30);
419 #if ROOT_VERSION_CODE < ROOT_VERSION(6,4,2)
420  h->SetBit(TH1::kCanRebin);
421 #else
422  h->SetCanExtend(TH1::kAllAxes);
423 #endif
424  std::array<std::string, 10> labels = {"PhysSel", "Evt Gen Name", "Trg class (acc)", "Trg class (rej)", "Cent", "vertex contr.", "Vz", "VzSPD", "SelPtHardBin", "MCOutlier"};
425  int i = 1;
426  for (auto label : labels) {
427  h->GetXaxis()->SetBinLabel(i, label.c_str());
428  i++;
429  }
430  h->GetYaxis()->SetTitle("counts");
431  fOutput->Add(h);
432  fHistograms["fHistEventRejection"] = h;
433 
434  h = new TH1F("fHistTriggerClasses","fHistTriggerClasses",3,0,3);
435 #if ROOT_VERSION_CODE < ROOT_VERSION(6,4,2)
436  h->SetBit(TH1::kCanRebin);
437 #else
438  h->SetCanExtend(TH1::kAllAxes);
439 #endif
440  fOutput->Add(h);
441  fHistograms["fHistTriggerClasses"] = h;
442 
443  h = new TH1F("fHistTriggerClassesNoSel","fHistTriggerClassesNoSel",3,0,3);
444 #if ROOT_VERSION_CODE < ROOT_VERSION(6,4,2)
445  h->SetBit(TH1::kCanRebin);
446 #else
447  h->SetCanExtend(TH1::kAllAxes);
448 #endif
449  fOutput->Add(h);
450  fHistograms["fHistTriggerClassesNoSel"] = h;
451 
452  h = new TH1F("fHistEventCount","fHistEventCount",2,0,2);
453  h->GetXaxis()->SetBinLabel(1,"Accepted");
454  h->GetXaxis()->SetBinLabel(2,"Rejected");
455  h->GetYaxis()->SetTitle("counts");
456  fOutput->Add(h);
457  fHistograms["fHistEventCount"] = h;
458 
459  // Finish setting up AliEventCuts
461  fAliEventCuts.AddQAplotsToList(fOutput);
462  }
463 
464  PostData(1, fOutput);
465 }
466 
468 {
469  if (eventSelected) {
470  if (fIsMonteCarlo) {
471  GetGeneralTH1("fHistEventsVsPtHard", true)->Fill(fPtHard);
472  GetGeneralTH1("fHistTrialsVsPtHard", true)->Fill(fPtHard, fNTrials);
473  GetGeneralTH1("fHistEventWeights", true)->Fill(fEventWeight);
474  GetGeneralTH2("fHistEventWeightsVsPtHard", true)->Fill(fPtHard, fEventWeight);
475  GetGeneralTH1("fHistXsectionDistribution", true)->Fill(fXsection);
476  GetGeneralTProfile("fHistXsection", true)->Fill(fPtHardBin, fXsection);
477  }
478 
479  GetGeneralTH1("fHistZVertex")->Fill(fVertex[2]);
480 
481  TH1* hCent = GetGeneralTH1("fHistCentrality");
482  if (hCent) hCent->Fill(fCent);
483 
484  TH1* hEventPlane = GetGeneralTH1("fHistEventPlane");
485  if (hEventPlane) hEventPlane->Fill(fEPV0);
486 
487  TH1* hTriggerClasses = GetGeneralTH1("fHistTriggerClasses");
488  for (auto fired_trg : fFiredTriggerClasses) hTriggerClasses->Fill(fired_trg.c_str(), 1);
489  }
490  else {
491  if (fIsMonteCarlo) {
492  GetGeneralTH1("fHistEventsVsPtHardNoSel", true)->Fill(fPtHard);
493  GetGeneralTH1("fHistTrialsVsPtHardNoSel", true)->Fill(fPtHard, fNTrials);
494  GetGeneralTH1("fHistEventWeightsNoSel", true)->Fill(fEventWeight);
495  GetGeneralTH2("fHistEventWeightsVsPtHardNoSel", true)->Fill(fPtHard, fEventWeight);
496  GetGeneralTH1("fHistXsectionDistributionNoSel", true)->Fill(fXsection);
497  GetGeneralTProfile("fHistXsectionNoSel", true)->Fill(fPtHardBin, fXsection);
498  }
499 
500  GetGeneralTH1("fHistZVertexNoSel", true)->Fill(fVertex[2]);
501 
502  TH1* hCent = GetGeneralTH1("fHistCentralityNoSel");
503  if (hCent) hCent->Fill(fCent);
504 
505  TH1* hEventPlane = GetGeneralTH1("fHistEventPlaneNoSel");
506  if (hEventPlane) hEventPlane->Fill(fEPV0);
507 
508  TH1* hTriggerClasses = GetGeneralTH1("fHistTriggerClassesNoSel", true);
509  for (auto fired_trg : fFiredTriggerClasses) hTriggerClasses->Fill(fired_trg.c_str(), 1);
510  }
511 
512  return kTRUE;
513 }
514 
516 {
517  if (fInhibit) {
518  AliWarningStream() << "The execution of this task is inhibited. Returning." << std::endl;
519  return;
520  }
521 
522  if (!fLocalInitialized) ExecOnce();
523 
524  if (!fLocalInitialized) return;
525 
526  if (!RetrieveEventObjects()) return;
527 
528  Bool_t eventSelected = IsEventSelected();
529 
531  if (eventSelected) {
532  GetGeneralTH1("fHistEventCount", true)->Fill("Accepted",1);
533  }
534  else {
535  GetGeneralTH1("fHistEventCount", true)->Fill("Rejected",1);
536  }
537 
538  FillGeneralHistograms(kFALSE);
539  if (eventSelected) FillGeneralHistograms(kTRUE);
540  }
541 
542  Bool_t runOk = kFALSE;
543  if (eventSelected || fEventSelectionAfterRun) runOk = Run();
544 
545  if (fCreateHisto && eventSelected && runOk) FillHistograms();
546 
547  if (fCreateHisto && fOutput) {
548  // information for this iteration of the UserExec in the container
549  PostData(1, fOutput);
550  }
551 }
552 
553 Bool_t AliAnalysisTaskEmcalLight::PythiaInfoFromFile(const char* currFile, Float_t &xsec, Float_t &trials, Int_t &pthard, Bool_t &useXsecFromHeader)
554 {
555 
556  TString file(currFile);
557  xsec = 0;
558  trials = 1;
559 
560  // Determine archive type
561  TString archivetype;
562  std::unique_ptr<TObjArray> walk(file.Tokenize("/"));
563  for(auto t : *walk){
564  TString &tok = static_cast<TObjString *>(t)->String();
565  if(tok.Contains(".zip")){
566  archivetype = tok;
567  Int_t pos = archivetype.Index(".zip");
568  archivetype.Replace(pos, archivetype.Length() - pos, "");
569  }
570  }
571  if(archivetype.Length()){
572  AliDebugStream(1) << "Auto-detected archive type " << archivetype << std::endl;
573  Ssiz_t pos1 = file.Index(archivetype,archivetype.Length(),0,TString::kExact);
574  Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
575  Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
576  file.Replace(pos+1,pos2-pos1,"");
577  } else {
578  // not an archive take the basename....
579  file.ReplaceAll(gSystem->BaseName(file.Data()),"");
580  }
581  AliDebugStream(1) << "File name: " << file << std::endl;
582 
583  // Build virtual file name
584  // Support for train tests
585  TString virtualFileName;
586  if(file.Contains("__alice")){
587  TString tmp(file);
588  Int_t pos = tmp.Index("__alice");
589  tmp.Replace(0, pos, "");
590  tmp.ReplaceAll("__", "/");
591  // cut out tag for archive and root file
592  // this needs a determin
593  std::unique_ptr<TObjArray> toks(tmp.Tokenize("/"));
594  TString tag = "_" + archivetype;
595  for(auto t : *toks){
596  TString &path = static_cast<TObjString *>(t)->String();
597  if(path.Contains(tag)){
598  Int_t posTag = path.Index(tag);
599  path.Replace(posTag, path.Length() - posTag, "");
600  }
601  virtualFileName += "/" + path;
602  }
603  } else {
604  virtualFileName = file;
605  }
606 
607  AliDebugStream(1) << "Physical file name " << file << ", virtual file name " << virtualFileName << std::endl;
608 
609  // Get the pt hard bin
610  TString strPthard(virtualFileName);
611 
612  /*
613  // Dead code - to be removed after testing phase
614  // Procedure will fail for everything else than the expected path name
615  strPthard.Remove(strPthard.Last('/'));
616  strPthard.Remove(strPthard.Last('/'));
617  if (strPthard.Contains("AOD")) strPthard.Remove(strPthard.Last('/'));
618  strPthard.Remove(0,strPthard.Last('/')+1);
619  if (strPthard.IsDec()) pthard = strPthard.Atoi();
620  else
621  AliWarningStream() << "Could not extract file number from path " << strPthard << std::endl;
622  */
623 
624  // New implementation : pattern matching
625  // Reason: Implementation valid only for old productions (new productions swap run number and pt-hard bin)
626  // Idea: Don't use the position in the string but the match different informations
627  // + Year clearly 2000+
628  // + Run number can be match to the one in the event
629  // + If we know it is not year or run number, it must be the pt-hard bin if we start from the beginning
630  // The procedure is only valid for the current implementations and unable to detect non-pt-hard bins
631  // It will also fail in case of arbitrary file names
632 
633  bool binfound = false;
634  std::unique_ptr<TObjArray> tokens(strPthard.Tokenize("/"));
635  for(auto t : *tokens) {
636  TString &tok = static_cast<TObjString *>(t)->String();
637  if(tok.IsDec()){
638  Int_t number = tok.Atoi();
639  if(number > 2000 && number < 3000){
640  // Year
641  continue;
642  } else if(number == fInputHandler->GetEvent()->GetRunNumber()){
643  // Run number
644  continue;
645  } else {
646  if(!binfound){
647  // the first number that is not one of the two must be the pt-hard bin
648  binfound = true;
649  pthard = number;
650  break;
651  }
652  }
653  }
654  }
655  if(!binfound) {
656  AliErrorStream() << "Could not extract file number from path " << strPthard << std::endl;
657  } else {
658  AliInfoStream() << "Auto-detecting pt-hard bin " << pthard << std::endl;
659  }
660 
661  AliInfoStream() << "File: " << file << std::endl;
662 
663  // problem that we cannot really test the existance of a file in a archive so we have to live with open error message from root
664  std::unique_ptr<TFile> fxsec(TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")));
665 
666  if (!fxsec) {
667  // next trial fetch the histgram file
668  fxsec = std::unique_ptr<TFile>(TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root")));
669  if (!fxsec){
670  AliErrorStream() << "Failed reading cross section from file " << file << std::endl;
671  useXsecFromHeader = true;
672  return kFALSE; // not a severe condition but inciate that we have no information
673  }
674  else {
675  // find the tlist we want to be independtent of the name so use the Tkey
676  TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
677  if (!key) return kFALSE;
678  TList *list = dynamic_cast<TList*>(key->ReadObj());
679  if (!list) return kFALSE;
680  TProfile *xSecHist = static_cast<TProfile*>(list->FindObject("h1Xsec"));
681  // check for failure
682  if(!xSecHist->GetEntries()) {
683  // No cross seciton information available - fall back to raw
684  AliErrorStream() << "No cross section information available in file " << fxsec->GetName() <<" - fall back to cross section in PYTHIA header" << std::endl;
685  useXsecFromHeader = true;
686  } else {
687  // Cross section histogram filled - take it from there
688  xsec = xSecHist->GetBinContent(1);
689  if(!xsec) AliErrorStream() << GetName() << ": Cross section 0 for file " << file << std::endl;
690  useXsecFromHeader = false;
691  }
692  trials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
693  }
694  } else { // no tree pyxsec.root
695  TTree *xtree = (TTree*)fxsec->Get("Xsection");
696  if (!xtree) return kFALSE;
697  UInt_t ntrials = 0;
698  Double_t xsection = 0;
699  xtree->SetBranchAddress("xsection",&xsection);
700  xtree->SetBranchAddress("ntrials",&ntrials);
701  xtree->GetEntry(0);
702  trials = ntrials;
703  xsec = xsection;
704  }
705  return kTRUE;
706 }
707 
709 {
710  if (!fIsPythia) return kTRUE;
711 
712  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
713  if (!tree) {
714  AliError(Form("%s - UserNotify: No current tree!",GetName()));
715  return kFALSE;
716  }
717 
718  Float_t xsection = 0;
719  Float_t trials = 0;
720  Int_t pthardbin = 0;
721  Bool_t useXsecFromHeader = false;
722 
723  TFile *curfile = tree->GetCurrentFile();
724  if (!curfile) {
725  AliError(Form("%s - UserNotify: No current file!",GetName()));
726  return kFALSE;
727  }
728 
729  TChain *chain = dynamic_cast<TChain*>(tree);
730  if (chain) tree = chain->GetTree();
731 
732  Int_t nevents = tree->GetEntriesFast();
733 
734  Bool_t res = PythiaInfoFromFile(curfile->GetName(), xsection, trials, pthardbin, useXsecFromHeader);
735 
736  fPtHardBin = pthardbin >= 0 ? pthardbin : 0;
737  fUseXsecFromHeader = useXsecFromHeader;
738 
739  if (!res) return kTRUE;
740 
742  GetGeneralTH1("fHistTrialsExternalFile", true)->Fill(fPtHardBin, trials);
743  if(!useXsecFromHeader) GetGeneralTProfile("fHistXsectionExternalFile", true)->Fill(fPtHardBin, xsection);
744  GetGeneralTH1("fHistEventsExternalFile", true)->Fill(fPtHardBin, nevents);
745  }
746 
747  return kTRUE;
748 }
749 
751 {
752  if (!InputEvent()) {
753  AliError(Form("%s: Could not retrieve event! Returning!", GetName()));
754  return;
755  }
756 
757  if (fNeedEmcalGeom && !fGeom) {
758  fGeom = AliEMCALGeometry::GetInstanceFromRunNumber(InputEvent()->GetRunNumber());
759  if (!fGeom) {
760  AliFatal(Form("%s: Can not get EMCal geometry instance. If you do not need the EMCal geometry, disable it by setting task->SetNeedEmcalGeometry(kFALSE).", GetName()));
761  return;
762  }
763  }
764 
765  if (fSwitchOffLHC15oFaultyBranches && dynamic_cast<AliAODEvent*>(InputEvent())) {
766  TTree *aodTree = AliAnalysisManager::GetAnalysisManager()->GetTree();
767  aodTree->SetBranchStatus("D0toKpi.fPx", 0);
768  aodTree->SetBranchStatus("D0toKpi.fPy", 0);
769  aodTree->SetBranchStatus("D0toKpi.fPz", 0);
770  aodTree->SetBranchStatus("D0toKpi.fd0", 0);
771  aodTree->SetBranchStatus("Charm3Prong.fPx", 0);
772  aodTree->SetBranchStatus("Charm3Prong.fPy", 0);
773  aodTree->SetBranchStatus("Charm3Prong.fPz", 0);
774  aodTree->SetBranchStatus("Charm3Prong.fd0", 0);
775  aodTree->SetBranchStatus("Dstar.fPx", 0);
776  aodTree->SetBranchStatus("Dstar.fPy", 0);
777  aodTree->SetBranchStatus("Dstar.fPz", 0);
778  aodTree->SetBranchStatus("Dstar.fd0", 0);
779  }
780 
781  //Load all requested track branches - each container knows name already
782  for (auto cont_it : fParticleCollArray) cont_it.second->SetArray(InputEvent());
783 
784  //Load all requested cluster branches - each container knows name already
785  for (auto cont_it : fClusterCollArray) cont_it.second->SetArray(InputEvent());
786 
787  if (!fCaloCellsName.IsNull() && !fCaloCells) {
788  fCaloCells = dynamic_cast<AliVCaloCells*>(InputEvent()->FindListObject(fCaloCellsName));
789  if (!fCaloCells) {
790  AliError(Form("%s: Could not retrieve cells %s!", GetName(), fCaloCellsName.Data()));
791  return;
792  }
793  }
794 
795  if (!fCaloTriggersName.IsNull() && !fCaloTriggers) {
796  fCaloTriggers = dynamic_cast<AliVCaloTrigger*>(InputEvent()->FindListObject(fCaloTriggersName));
797  if (!fCaloTriggers) {
798  AliError(Form("%s: Could not retrieve calo triggers %s!", GetName(), fCaloTriggersName.Data()));
799  return;
800  }
801  }
802 
803  if (!fCaloTriggerPatchInfoName.IsNull() && !fTriggerPatchInfo) {
804  fTriggerPatchInfo = GetArrayFromEvent(fCaloTriggerPatchInfoName.Data(),"AliEMCALTriggerPatchInfo");
805  if (!fTriggerPatchInfo) {
806  AliError(Form("%s: Could not retrieve calo trigger patch info %s!", GetName(), fCaloTriggerPatchInfoName.Data()));
807  return;
808  }
809 
810  }
811 
812  fLocalInitialized = kTRUE;
813 }
814 
816 {
817  if (fForceBeamType != kNA)
818  return fForceBeamType;
819 
820  AliESDEvent *esd = dynamic_cast<AliESDEvent*>(InputEvent());
821  if (esd) {
822  const AliESDRun *run = esd->GetESDRun();
823  TString beamType = run->GetBeamType();
824  if (beamType == "p-p")
825  return kpp;
826  else if (beamType == "A-A")
827  return kAA;
828  else if (beamType == "p-A")
829  return kpA;
830  else
831  return kNA;
832  } else {
833  Int_t runNumber = InputEvent()->GetRunNumber();
834  // All run number ranges taken from the RCT
835  if ((runNumber >= 136833 && runNumber <= 139517) || // LHC10h
836  (runNumber >= 167693 && runNumber <= 170593) || // LHC11h
837  (runNumber >= 244824 && runNumber <= 246994)) { // LHC15o
838  return kAA;
839  } else if ((runNumber >= 188356 && runNumber <= 188366) || // LHC12g
840  (runNumber >= 195164 && runNumber <= 197388) || // LHC13b-f
841  (runNumber >= 265015 && runNumber <= 267166)) { // LHC16q-t
842  return kpA;
843  } else {
844  return kpp;
845  }
846  }
847 }
848 
850  if(!IsTriggerSelected()) return false;
852  if(!CheckMCOutliers()) return false;
853  return fAliEventCuts.AcceptEvent(fInputEvent);
854 }
855 
857 {
858  TH1* hEventRejection = GetGeneralTH1("fHistEventRejection", true);
859 
861  if (fGeneralHistograms) hEventRejection->Fill("PhysSel",1);
862  return kFALSE;
863  }
864 
865  if (!fSelectGeneratorName.IsNull() && !fGeneratorName.IsNull()) {
866  if (!fGeneratorName.Contains(fSelectGeneratorName)) {
867  if (fGeneralHistograms) hEventRejection->Fill("Evt Gen Name",1);
868  return kFALSE;
869  }
870  }
871 
872  if (fMinCent < fMaxCent && fMaxCent > 0) {
873  if (fCent < fMinCent || fCent > fMaxCent) {
874  if (fGeneralHistograms) hEventRejection->Fill("Cent",1);
875  return kFALSE;
876  }
877  }
878 
879  if (fNVertCont < fMinNVertCont) {
880  if (fGeneralHistograms) hEventRejection->Fill("vertex contr.",1);
881  return kFALSE;
882  }
883 
884  if (fMinVz < fMaxVz) {
885  if (fVertex[2] < fMinVz || fVertex[2] > fMaxVz) {
886  if (fGeneralHistograms) hEventRejection->Fill("Vz",1);
887  return kFALSE;
888  }
889  }
890 
891  if (fMaxVzDiff >= 0) {
892  if (fNVertSPDCont > 0) {
893  Double_t vzSPD = fVertexSPD[2];
894  Double_t dvertex = TMath::Abs(fVertex[2] - vzSPD);
895  //if difference larger than fZvertexDiff
896  if (dvertex > fMaxVzDiff) {
897  if (fGeneralHistograms) hEventRejection->Fill("VzSPD",1);
898  return kFALSE;
899  }
900  }
901  }
902 
903  if (fMinPtHard >= 0 && fPtHard < fMinPtHard) {
904  if (fGeneralHistograms) hEventRejection->Fill("SelPtHardBin",1);
905  return kFALSE;
906  }
907 
908  if (fMaxPtHard >= 0 && fPtHard >= fMaxPtHard) {
909  if (fGeneralHistograms) hEventRejection->Fill("SelPtHardBin",1);
910  return kFALSE;
911  }
912 
914  if (fGeneralHistograms) hEventRejection->Fill("SelPtHardBin",1);
915  return kFALSE;
916  }
917 
918  // Reject filter for MC data
919  if (!CheckMCOutliers()) {
920  if (fGeneralHistograms) hEventRejection->Fill("MCOutlier",1);
921  return kFALSE;
922  }
923 
924  return kTRUE;
925 }
926 
928  TH1* hEventRejection = GetGeneralTH1("fHistEventRejection", true);
929  Bool_t acceptedTrgClassFound = kFALSE;
930  if (fAcceptedTriggerClasses.size() > 0) {
931  for (const auto &acc_trg : fAcceptedTriggerClasses) {
932  std::string teststring(acc_trg);
933  bool fullmatch(false);
934  auto posexact = acc_trg.find("EXACT");
935  if(posexact != std::string::npos) {
936  fullmatch = true;
937  teststring.erase(posexact, 5);
938  }
939  for (const auto &fired_trg : fFiredTriggerClasses) {
940  bool classmatch = fullmatch ? teststring == fired_trg : fired_trg.find(teststring) != std::string::npos;
941  if (classmatch) {
942  acceptedTrgClassFound = kTRUE;
943  break;
944  }
945  }
946  if (acceptedTrgClassFound) break;
947  }
948 
949  if (!acceptedTrgClassFound) {
950  if (fGeneralHistograms) hEventRejection->Fill("Trg class (acc)",1);
951  return kFALSE;
952  }
953  }
954 
955  if (fRejectedTriggerClasses.size() > 0) {
956  for (const auto &rej_trg : fRejectedTriggerClasses) {
957  std::string teststring(rej_trg);
958  bool fullmatch(false);
959  auto posexact = rej_trg.find("EXACT");
960  if(posexact != std::string::npos) {
961  fullmatch = true;
962  teststring.erase(posexact, 5);
963  }
964  for (const auto &fired_trg : fFiredTriggerClasses) {
965  bool classmatch = fullmatch ? teststring == fired_trg : fired_trg.find(teststring) != std::string::npos;
966  if (classmatch) {
967  if (fGeneralHistograms) hEventRejection->Fill("Trg class (rej)",1);
968  return kFALSE;
969  }
970  }
971  }
972  }
973  return kTRUE;
974 }
975 
976 TClonesArray *AliAnalysisTaskEmcalLight::GetArrayFromEvent(const char *name, const char *clname)
977 {
978  TClonesArray *arr = 0;
979  TString sname(name);
980  if (!sname.IsNull()) {
981  arr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(sname));
982  if (!arr) {
983  AliWarning(Form("%s: Could not retrieve array with name %s!", GetName(), name));
984  return 0;
985  }
986  } else {
987  return 0;
988  }
989 
990  if (!clname)
991  return arr;
992 
993  TString objname(arr->GetClass()->GetName());
994  TClass cls(objname);
995  if (!cls.InheritsFrom(clname)) {
996  AliWarning(Form("%s: Objects of type %s in %s are not inherited from %s!",
997  GetName(), cls.GetName(), name, clname));
998  return 0;
999  }
1000  return arr;
1001 }
1002 
1004 {
1005  fVertex[0] = 0;
1006  fVertex[1] = 0;
1007  fVertex[2] = 0;
1008  fNVertCont = 0;
1009 
1010  fVertexSPD[0] = 0;
1011  fVertexSPD[1] = 0;
1012  fVertexSPD[2] = 0;
1013  fNVertSPDCont = 0;
1014 
1015  fFiredTriggerClasses.clear();
1016  std::stringstream firedClasses(InputEvent()->GetFiredTriggerClasses().Data());
1017  while (firedClasses.good()) {
1018  std::string trgClass;
1019  firedClasses >> trgClass;
1020  if (!trgClass.empty()) fFiredTriggerClasses.push_back(trgClass);
1021  }
1022 
1023  if (fDataType == kESD) {
1024  fFiredTriggerBitMap = static_cast<AliInputEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler())->IsEventSelected();
1025  }
1026  else {
1027  fFiredTriggerBitMap = static_cast<AliVAODHeader*>(InputEvent()->GetHeader())->GetOfflineTrigger();
1028  }
1029 
1030  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
1031  if (vert) {
1032  vert->GetXYZ(fVertex);
1033  fNVertCont = vert->GetNContributors();
1034  }
1035 
1036  const AliVVertex *vertSPD = InputEvent()->GetPrimaryVertexSPD();
1037  if (vertSPD) {
1038  vertSPD->GetXYZ(fVertexSPD);
1039  fNVertSPDCont = vertSPD->GetNContributors();
1040  }
1041 
1042  fBeamType = GetBeamType();
1043 
1044  fCent = 99;
1045  fCentBin = -1;
1046  fEPV0 = -999;
1047  fEPV0A = -999;
1048  fEPV0C = -999;
1049 
1051  // New centrality estimation (AliMultSelection)
1052  // See https://twiki.cern.ch/twiki/bin/viewauth/ALICE/AliMultSelectionCalibStatus for calibration status period-by-period)
1053  AliMultSelection *MultSelection = static_cast<AliMultSelection*>(InputEvent()->FindListObject("MultSelection"));
1054  if (MultSelection) {
1055  fCent = MultSelection->GetMultiplicityPercentile(fCentEst.Data());
1056  }
1057  else {
1058  if(fWarnMissingCentrality) AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
1059  }
1060  }
1061  else if (fCentralityEstimation == kOldCentrality) {
1062  // Old centrality estimation (AliCentrality, works only on Run-1 PbPb and pPb)
1063  AliCentrality *aliCent = InputEvent()->GetCentrality();
1064  if (aliCent) {
1065  fCent = aliCent->GetCentralityPercentile(fCentEst.Data());
1066  }
1067  else {
1068  if(fWarnMissingCentrality) AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
1069  }
1070  }
1071  if (!fCentBins.empty() && fCentralityEstimation != kNoCentrality) {
1072  for (auto cent_it = fCentBins.begin(); cent_it != fCentBins.end() - 1; cent_it++) {
1073  if (fCent >= *cent_it && fCent < *(cent_it+1)) fCentBin = cent_it - fCentBins.begin();
1074  }
1075  }
1076  else {
1077  fCentBin = 0;
1078  }
1079 
1080  if (fBeamType == kAA || fBeamType == kpA ) {
1081  AliEventplane *aliEP = InputEvent()->GetEventplane();
1082  if (aliEP) {
1083  fEPV0 = aliEP->GetEventplane("V0" ,InputEvent());
1084  fEPV0A = aliEP->GetEventplane("V0A",InputEvent());
1085  fEPV0C = aliEP->GetEventplane("V0C",InputEvent());
1086  } else {
1087  AliWarning(Form("%s: Could not retrieve event plane information!", GetName()));
1088  }
1089  }
1090 
1091  if (fIsMonteCarlo && MCEvent()) {
1092  AliGenEventHeader* header = MCEvent()->GenEventHeader();
1093  if (fMCEventHeaderName.IsNull()) {
1094  fMCHeader = header;
1095  }
1096  else {
1097  if (header->InheritsFrom(fMCEventHeaderName)) {
1098  fMCHeader = header;
1099  }
1100  else if (header->InheritsFrom("AliGenCocktailEventHeader")) {
1101  AliGenCocktailEventHeader* cocktailHeader = static_cast<AliGenCocktailEventHeader*>(header);
1102  TList* headers = cocktailHeader->GetHeaders();
1103  for (auto obj : *headers) { // @suppress("Symbol is not resolved")
1104  if (obj->InheritsFrom(fMCEventHeaderName)){
1105  fMCHeader = static_cast<AliGenEventHeader*>(obj);
1106  break;
1107  }
1108  }
1109  }
1110  }
1111  if (fMCHeader) {
1112  fEventWeight = fMCHeader->EventWeight();
1113  if (fIsPythia) {
1114  fPythiaHeader = static_cast<AliGenPythiaEventHeader*>(fMCHeader);
1115  fPtHard = fPythiaHeader->GetPtHard();
1116  fXsection = fPythiaHeader->GetXsection();
1117  fNTrials = fPythiaHeader->Trials();
1118  if(fUseXsecFromHeader) GetGeneralTProfile("fHistXsectionExternalFile", true)->Fill(fPtHardBin, fXsection);
1119  }
1120  }
1121  }
1122 
1123  for (auto cont_it : fParticleCollArray) cont_it.second->NextEvent(InputEvent());
1124  for (auto cont_it : fClusterCollArray) cont_it.second->NextEvent(InputEvent());
1125 
1126  return kTRUE;
1127 }
1128 
1130 {
1131  if (branchName.size() == 0) return 0;
1132 
1133  AliParticleContainer* cont = 0;
1134 
1135 #if ROOT_VERSION_CODE > ROOT_VERSION(6,10,0)
1136 #define EMCAL_STRINGVIEW_NONCONST std::string_view
1137 #else
1138 #define EMCAL_STRINGVIEW_NONCONST std::string
1139 #endif
1140 
1141  if (branchName == EMCAL_STRINGVIEW_NONCONST("tracks") || branchName == EMCAL_STRINGVIEW_NONCONST("Tracks")) cont = new AliTrackContainer(branchName.data());
1142  else if (branchName == EMCAL_STRINGVIEW_NONCONST("mcparticles")) cont = new AliMCParticleContainer(branchName.data());
1143  else cont = new AliParticleContainer(branchName.data());
1144 
1145  if (contName.size() > 0) cont->SetName(contName.data());
1146 
1147  AdoptParticleContainer(cont);
1148 
1149  return cont;
1150 }
1151 
1153 {
1154  if (branchName.size() == 0) return 0;
1155 
1156  AliClusterContainer* cont = new AliClusterContainer(branchName.data());
1157 
1158  if (contName.size() > 0) cont->SetName(contName.data());
1159 
1160  AdoptClusterContainer(cont);
1161 
1162  return cont;
1163 }
1164 
1166 {
1167  std::map<std::string, AliParticleContainer*>::const_iterator cont_it = fParticleCollArray.find(std::string(name));
1168  if (cont_it != fParticleCollArray.end()) return cont_it->second;
1169  else return nullptr;
1170 }
1171 
1173 {
1174  std::map<std::string, AliClusterContainer*>::const_iterator cont_it = fClusterCollArray.find(std::string(name));
1175  if (cont_it != fClusterCollArray.end()) return cont_it->second;
1176  else return nullptr;
1177 }
1178 
1180 {
1181  if (!(InputEvent()->FindListObject(obj->GetName()))) {
1182  InputEvent()->AddObject(obj);
1183  }
1184  else {
1185  if (!attempt) {
1186  AliFatal(Form("%s: Container with name %s already present. Aborting", GetName(), obj->GetName()));
1187  }
1188  }
1189 }
1190 
1192 {
1193 
1194  if (!fGeom) {
1195  AliWarning(Form("%s - AliAnalysisTaskEmcalBase::IsTrackInEmcalAcceptance - Geometry is not available!", GetName()));
1196  return kFALSE;
1197  }
1198 
1199  Double_t minPhi = fGeom->GetArm1PhiMin() - edges;
1200  Double_t maxPhi = fGeom->GetArm1PhiMax() + edges;
1201 
1202  if (part->Phi() > minPhi && part->Phi() < maxPhi) {
1203  return kTRUE;
1204  }
1205  else {
1206  return kFALSE;
1207  }
1208 }
1209 
1211 {
1212  axis->SetBinLabel(1, "NullObject");
1213  axis->SetBinLabel(2, "Pt");
1214  axis->SetBinLabel(3, "Acceptance");
1215  axis->SetBinLabel(4, "MCLabel");
1216  axis->SetBinLabel(5, "BitMap");
1217  axis->SetBinLabel(6, "HF cut");
1218  axis->SetBinLabel(7, "Bit6");
1219  axis->SetBinLabel(8, "NotHybridTrack");
1220  axis->SetBinLabel(9, "MCFlag");
1221  axis->SetBinLabel(10, "MCGenerator");
1222  axis->SetBinLabel(11, "ChargeCut");
1223  axis->SetBinLabel(12, "MinDistanceTPCSectorEdge");
1224  axis->SetBinLabel(13, "Bit12");
1225  axis->SetBinLabel(14, "IsEMCal");
1226  axis->SetBinLabel(15, "Time");
1227  axis->SetBinLabel(16, "Energy");
1228  axis->SetBinLabel(17, "ExoticCut");
1229  axis->SetBinLabel(18, "Bit17");
1230  axis->SetBinLabel(19, "Area");
1231  axis->SetBinLabel(20, "AreaEmc");
1232  axis->SetBinLabel(21, "ZLeadingCh");
1233  axis->SetBinLabel(22, "ZLeadingEmc");
1234  axis->SetBinLabel(23, "NEF");
1235  axis->SetBinLabel(24, "MinLeadPt");
1236  axis->SetBinLabel(25, "MaxTrackPt");
1237  axis->SetBinLabel(26, "MaxClusterPt");
1238  axis->SetBinLabel(27, "Flavour");
1239  axis->SetBinLabel(28, "TagStatus");
1240  axis->SetBinLabel(29, "MinNConstituents");
1241  axis->SetBinLabel(30, "Bit29");
1242  axis->SetBinLabel(31, "Bit30");
1243  axis->SetBinLabel(32, "Bit31");
1244 }
1245 
1246 Double_t AliAnalysisTaskEmcalLight::GetParallelFraction(AliVParticle* part1, AliVParticle* part2)
1247 {
1248  TVector3 vect1(part1->Px(), part1->Py(), part1->Pz());
1249  TVector3 vect2(part2->Px(), part2->Py(), part2->Pz());
1250  Double_t z = (vect1 * vect2) / (vect2 * vect2);
1251  return z;
1252 }
1253 
1254 Double_t AliAnalysisTaskEmcalLight::GetParallelFraction(const TVector3& vect1, AliVParticle* part2)
1255 {
1256  TVector3 vect2(part2->Px(), part2->Py(), part2->Pz());
1257  Double_t z = (vect1 * vect2) / (vect2 * vect2);
1258  return z;
1259 }
1260 
1261 void AliAnalysisTaskEmcalLight::GetEtaPhiDiff(const AliVTrack *t, const AliVCluster *v, Double_t &phidiff, Double_t &etadiff)
1262 {
1263  phidiff = 999;
1264  etadiff = 999;
1265 
1266  if (!t||!v) return;
1267 
1268  Double_t veta = t->GetTrackEtaOnEMCal();
1269  Double_t vphi = t->GetTrackPhiOnEMCal();
1270 
1271  Float_t pos[3] = {0};
1272  v->GetPosition(pos);
1273  TVector3 cpos(pos);
1274  Double_t ceta = cpos.Eta();
1275  Double_t cphi = cpos.Phi();
1276  etadiff=veta-ceta;
1277  phidiff=TVector2::Phi_mpi_pi(vphi-cphi);
1278 }
1279 
1281 {
1282  Byte_t ret = 0;
1283  if (t->TestBit(BIT(22)) && !t->TestBit(BIT(23)))
1284  ret = 1;
1285  else if (!t->TestBit(BIT(22)) && t->TestBit(BIT(23)))
1286  ret = 2;
1287  else if (t->TestBit(BIT(22)) && t->TestBit(BIT(23)))
1288  ret = 3;
1289  return ret;
1290 }
1291 
1292 Byte_t AliAnalysisTaskEmcalLight::GetTrackType(const AliAODTrack *aodTrack, UInt_t filterBit1, UInt_t filterBit2)
1293 {
1294 
1295  Int_t res = 0;
1296 
1297  if (aodTrack->TestFilterBit(filterBit1)) {
1298  res = 0;
1299  }
1300  else if (aodTrack->TestFilterBit(filterBit2)) {
1301  if ((aodTrack->GetStatus()&AliVTrack::kITSrefit)!=0) {
1302  res = 1;
1303  }
1304  else {
1305  res = 2;
1306  }
1307  }
1308  else {
1309  res = 3;
1310  }
1311 
1312  return res;
1313 }
1314 
1316 {
1317  EBeamType_t b = kpp;
1318  if ((runnumber >= 136833 && runnumber <= 139517) || // LHC10h Run-1 (Pb-Pb)
1319  (runnumber >= 167693 && runnumber <= 170593) || // LHC11h Run-1 (Pb-Pb)
1320  (runnumber >= 244824 && runnumber <= 246994) //|| // LHC15o Run-2 (Pb-Pb)
1321  //(runnumber >= 295581 && runnumber <= 297624) // LHC18q+r Run-2 (Pb-Pb)
1322  )
1323  {
1324  b = kAA;
1325  }
1326  else if ((runnumber > 188356 && runnumber <= 188503) || // LHC12g Run-1 (p-Pb pilot)
1327  (runnumber >= 195164 && runnumber <= 197388) || // LHC13b,c,d,e,f Run-1 (p-Pb)
1328  (runnumber >= 265077 && runnumber <= 267166)) { // LHC16 Run-2 (p-Pb)
1329  b = kpA;
1330  }
1331  return b;
1332 }
1333 
1335 {
1336  if (!fPythiaHeader || !fMCRejectFilter) return kTRUE;
1337 
1338  // Condition 1: Pythia jet / pT-hard > factor
1339  if (fPtHardAndJetPtFactor > 0.) {
1340  AliTLorentzVector jet;
1341 
1342  Int_t nTriggerJets = fPythiaHeader->NTriggerJets();
1343 
1344  AliDebug(1,Form("Njets: %d, pT Hard %f",nTriggerJets, fPtHard));
1345 
1346  Float_t tmpjet[]={0,0,0,0};
1347  for (Int_t ijet = 0; ijet< nTriggerJets; ijet++) {
1348  fPythiaHeader->TriggerJet(ijet, tmpjet);
1349 
1350  jet.SetPxPyPzE(tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3]);
1351 
1352  AliDebug(1,Form("jet %d; pycell jet pT %f",ijet, jet.Pt()));
1353 
1354  //Compare jet pT and pt Hard
1355  if (jet.Pt() > fPtHardAndJetPtFactor * fPtHard) {
1356  AliInfo(Form("Reject jet event with : pT Hard %2.2f, pycell jet pT %2.2f, rejection factor %1.1f\n", fPtHard, jet.Pt(), fPtHardAndJetPtFactor));
1357  return kFALSE;
1358  }
1359  }
1360  }
1361  // end condition 1
1362 
1363  // Condition 2 : Reconstructed EMCal cluster pT / pT-hard > factor
1364  if (fPtHardAndClusterPtFactor > 0.) {
1365  AliClusterContainer* mccluscont = fClusterCollArray.begin()->second;
1366  if ((Bool_t)mccluscont) {
1367  for (auto cluster : mccluscont->all()) {// Not cuts applied ; use accept for cuts
1368  Float_t ecluster = cluster->E();
1369 
1370  if (ecluster > (fPtHardAndClusterPtFactor * fPtHard)) {
1371  AliInfo(Form("Reject : ecluster %2.2f, calo %d, factor %2.2f, ptHard %f",ecluster,cluster->GetType(),fPtHardAndClusterPtFactor,fPtHard));
1372  return kFALSE;
1373  }
1374  }
1375  }
1376  }
1377  // end condition 2
1378 
1379  // condition 3 : Reconstructed track pT / pT-hard >factor
1380  std::vector<AliMCParticleContainer *> mcpcont;
1381  for(auto cont : fParticleCollArray) {
1382  AliMCParticleContainer *mccont = dynamic_cast<AliMCParticleContainer *>(cont.second);
1383  if(mccont) mcpcont.push_back(mccont);
1384  }
1385  if (fPtHardAndTrackPtFactor > 0.) {
1386  AliMCParticleContainer* mcpartcont = *mcpcont.begin();
1387  if ((Bool_t)mcpartcont) {
1388  for (auto mctrack : mcpartcont->all()) {// Not cuts applied ; use accept for cuts
1389  Float_t trackpt = mctrack->Pt();
1390  if (trackpt > (fPtHardAndTrackPtFactor * fPtHard) ) {
1391  AliInfo(Form("Reject : track %2.2f, factor %2.2f, ptHard %f", trackpt, fPtHardAndTrackPtFactor, fPtHard));
1392  return kFALSE;
1393  }
1394  }
1395  }
1396  }
1397  // end condition 3
1398 
1399  return kTRUE;
1400 }
1401 
1403 {
1404  Double_t dphi = -999;
1405  const Double_t tpi = TMath::TwoPi();
1406 
1407  if (phia < 0) phia += tpi;
1408  else if (phia > tpi) phia -= tpi;
1409  if (phib < 0) phib += tpi;
1410  else if (phib > tpi) phib -= tpi;
1411  dphi = phib - phia;
1412  if (dphi < rangeMin) dphi += tpi;
1413  else if (dphi > rangeMax) dphi -= tpi;
1414 
1415  return dphi;
1416 }
1417 
1418 void AliAnalysisTaskEmcalLight::GenerateFixedBinArray(int n, double min, double max, std::vector<double>& array, bool last)
1419 {
1420  double binWidth = (max - min) / n;
1421  double v = min;
1422  if (last) n++;
1423  for (int i = 0; i < n; i++) {
1424  array.push_back(v);
1425  v += binWidth;
1426  }
1427 }
1428 
1429 std::vector<double> AliAnalysisTaskEmcalLight::GenerateFixedBinArray(int n, double min, double max, bool last)
1430 {
1431  std::vector<double> array;
1432  GenerateFixedBinArray(n, min, max, array, last);
1433  return array;
1434 }
1435 
1436 void AliAnalysisTaskEmcalLight::GenerateLogFixedBinArray(int n, double min, double max, std::vector<double>& array, bool last)
1437 {
1438  if (min <= 0 || max < min) {
1439  AliErrorClassStream() << "Cannot generate a log scale fixed-bin array with limits " << min << ", " << max << std::endl;
1440  return;
1441  }
1442  double binWidth = std::pow(max / min, 1.0 / n);
1443  double v = min;
1444  if (last) n++;
1445  for (int i = 0; i < n; i++) {
1446  array.push_back(v);
1447  v *= binWidth;
1448  }
1449 }
1450 
1451 std::vector<double> AliAnalysisTaskEmcalLight::GenerateLogFixedBinArray(int n, double min, double max, bool last)
1452 {
1453  std::vector<double> array;
1454  GenerateLogFixedBinArray(n, min, max, array, last);
1455  return array;
1456 }
1457 
1458 
1459 TH1* AliAnalysisTaskEmcalLight::GetGeneralTH1(const char* name, bool warn)
1460 {
1461  auto search = fHistograms.find(name);
1462  if (search != fHistograms.end()) {
1463  return search->second;
1464  }
1465  else {
1466  if (warn) AliErrorStream() << "Could not find histogram '" << name << "'" << std::endl;
1467  return nullptr;
1468  }
1469 }
1470 
1471 TH2* AliAnalysisTaskEmcalLight::GetGeneralTH2(const char* name, bool warn)
1472 {
1473  return static_cast<TH2*>(GetGeneralTH1(name, warn));
1474 }
1475 
1476 TProfile* AliAnalysisTaskEmcalLight::GetGeneralTProfile(const char* name, bool warn)
1477 {
1478  return static_cast<TProfile*>(GetGeneralTH1(name, warn));
1479 }
1480 
1482 {
1483  fIsPythia = i;
1484  if (fIsPythia) {
1485  fIsMonteCarlo = kTRUE;
1486  fMCEventHeaderName = "AliGenPythiaEventHeader";
1487  }
1488  else {
1489  if (fMCEventHeaderName == "AliGenPythiaEventHeader") {
1490  fMCEventHeaderName = "";
1491  }
1492  }
1493 }
1494 
1496 {
1497  TClass gen_header_class(name);
1498  if (gen_header_class.InheritsFrom("AliGenEventHeader")) {
1499  fMCEventHeaderName = name;
1500  }
1501  else {
1502  AliWarningStream() << "Class name '" << name << "' does not inherit from 'AliGenEventHeader'. Not setting it." << std::endl;
1503  }
1504 }
Bool_t fSwitchOffLHC15oFaultyBranches
Switch off faulty tree branches in LHC15o AOD trees.
Float_t fPtHardAndJetPtFactor
Factor between ptHard and jet pT to reject/accept event.
Double_t fVertexSPD[3]
!event Svertex
TString fCaloTriggersName
name of calo triggers collection
AliEMCALGeometry * fGeom
!emcal geometry
EBeamType_t fBeamType
!event beam type
Bool_t fInhibit
!inhibit execution of the task
double Double_t
Definition: External.C:58
TClonesArray * GetArrayFromEvent(const char *name, const char *clname=0)
Definition: External.C:236
void SetName(const char *n)
Set the name of the class of the objets inside the underlying array.
Float_t fXsection
!x-section from pythia header
#define EMCAL_STRINGVIEW_NONCONST
Double_t fEPV0A
!event plane V0A
TString fSelectGeneratorName
Selects only events produced by a generator that has a name containing a string.
Double_t fEPV0
!event plane V0
Bool_t fIsMonteCarlo
if it is a MC production
Container with name, TClonesArray and cuts for particles.
virtual Bool_t IsTriggerSelected()
Perform trigger selection.
TSystem * gSystem
Bool_t fMCRejectFilter
enable the filtering of events by tail rejection
TString fCaloTriggerPatchInfoName
trigger patch info array name
Double_t fMinVz
min vertex for event selection
New centrality estimation (AliMultSelection, see https://twiki.cern.ch/twiki/bin/viewauth/ALICE/AliMu...
AliClusterContainer * AddClusterContainer(EMCAL_STRINGVIEW branchName, EMCAL_STRINGVIEW contName="")
TString fCentEst
name of the centrality estimator
UInt_t fTriggerSelectionBitMap
trigger selection bit map
std::vector< double > fCentBins
how many centrality bins
Bool_t fEventSelectionAfterRun
If kTRUE, the event selection is performed after Run() but before FillHistograms() ...
static Byte_t GetTrackType(const AliVTrack *t)
static std::vector< double > GenerateFixedBinArray(int n, double min, double max, bool last=true)
TString fCaloCellsName
name of calo cell collection
Bool_t IsTrackInEmcalAcceptance(AliVParticle *part, Double_t edges=0.9) const
void SetMCEventHeaderName(const char *name)
std::set< std::string > fRejectedTriggerClasses
list of accepted trigger classes
TH1 * GetGeneralTH1(const char *name, bool warn=false)
Double_t fMaxMinimumBiasPtHard
maximum pt hard for the minimum bias pt hard bin (MC)
Container for particles within the EMCAL framework.
static void GetEtaPhiDiff(const AliVTrack *t, const AliVCluster *v, Double_t &phidiff, Double_t &etadiff)
void AddObjectToEvent(TObject *obj, Bool_t attempt=kFALSE)
AliParticleContainer * AddParticleContainer(EMCAL_STRINGVIEW branchName, EMCAL_STRINGVIEW contName="")
std::map< std::string, TH1 * > fHistograms
!general QA histograms
static Double_t fgkEMCalDCalPhiDivide
phi value used to distinguish between DCal and EMCal
EBeamType_t fForceBeamType
forced beam type
std::vector< std::string > fFiredTriggerClasses
!trigger classes fired by the current event
Bool_t fCreateHisto
whether or not create histograms
Bool_t PythiaInfoFromFile(const char *currFile, Float_t &fXsec, Float_t &fTrials, Int_t &pthard, Bool_t &useXsecFromHeader)
int Int_t
Definition: External.C:63
TH2 * GetGeneralTH2(const char *name, bool warn=false)
Double_t fMinPtHard
select minimum pt hard (MC)
unsigned int UInt_t
Definition: External.C:33
TString fMCEventHeaderName
Looks for MC event properties in a particular MC event type (useful for a MC cocktail production) ...
float Float_t
Definition: External.C:68
Base task in the EMCAL framework (lighter version of AliAnalysisTaskEmcal)
Float_t fPtHardAndClusterPtFactor
Factor between ptHard and cluster pT to reject/accept event.
Double_t fMaxVzDiff
upper limit for distance between primary and SPD vertex
Bool_t fUsePtHardBinScaling
Apply pt-hard bin scaling (in case AliEmcalList is used to handle the output)
Double_t fVertex[3]
!event vertex
Double_t fMaximumEventWeight
Minimum event weight for the related bookkeping histogram.
ECentralityEstimation_t fCentralityEstimation
Centrality estimation.
void SetUseScaling(Bool_t val)
Definition: AliEmcalList.h:31
Int_t fPtHardBin
!event pt hard bin
void AdoptClusterContainer(AliClusterContainer *cont)
Bool_t fIsPythia
if it is a PYTHIA production
Int_t fCentBin
!event centrality bin
Double_t fMaxCent
max centrality for event selection
ULong_t fFiredTriggerBitMap
!bit map of fired triggers
Double_t fMinNVertCont
minumum number of vertex contributors
Int_t fNVertSPDCont
!event SPD vertex number of contributors
std::set< std::string > fAcceptedTriggerClasses
list of accepted trigger classes
static Double_t DeltaPhi(Double_t phia, Double_t phib, Double_t rMin=-TMath::Pi()/2, Double_t rMax=3 *TMath::Pi()/2)
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
Double_t fMinimumEventWeight
Minimum event weight for the related bookkeping histogram.
Enhanced TList-derived class that implements correct merging for pt_hard binned production.
Definition: AliEmcalList.h:25
Bool_t fLocalInitialized
!whether or not the task has been already initialized
AliGenEventHeader * fMCHeader
!event MC header
static EBeamType_t BeamTypeFromRunNumber(Int_t runnumber)
Bool_t fWarnMissingCentrality
!switch for verbosity in centrality information
Definition: External.C:220
TProfile * GetGeneralTProfile(const char *name, bool warn=false)
Bool_t fUseAliEmcalList
Use AliEmcalList as output object.
std::map< std::string, AliParticleContainer * > fParticleCollArray
particle/track collection array
TFile * file
TList with histograms for a given trigger.
Int_t nevents[nsamples]
#define EMCAL_STRINGVIEW
Int_t fNVertCont
!event vertex number of contributors
static std::vector< double > GenerateLogFixedBinArray(int n, double min, double max, bool last=true)
Bool_t fGeneralHistograms
whether or not it should fill some general histograms
AliClusterContainer * GetClusterContainer(EMCAL_STRINGVIEW name) const
Int_t GetRunNumber(TString)
Definition: PlotMuonQA.C:2235
AliGenPythiaEventHeader * fPythiaHeader
!event Pythia header
const char Option_t
Definition: External.C:48
Float_t fPtHardAndTrackPtFactor
Factor between ptHard and track pT to reject/accept event.
Bool_t fUseXsecFromHeader
!Switch for using cross section from header (if not found in pythia file)
AliEventCuts fAliEventCuts
Event cut object.
AliAnalysisTaskEmcalLight()
Default constructor.
bool Bool_t
Definition: External.C:53
TTree * tree
std::map< std::string, AliClusterContainer * > fClusterCollArray
cluster collection array
EBeamType_t
Switch for the beam type.
AliVCaloTrigger * fCaloTriggers
!calo triggers
Container structure for EMCAL clusters.
Container for MC-true particles within the EMCAL framework.
TString fGeneratorName
!name of the MC generator used to produce the current event (only AOD)
void AdoptParticleContainer(AliParticleContainer *cont)
virtual Bool_t FillGeneralHistograms(Bool_t eventSelected)
Double_t fMaxVz
max vertex for event selection
Double_t fCent
!event centrality
AliParticleContainer * GetParticleContainer(EMCAL_STRINGVIEW name) const
Double_t fEPV0C
!event plane V0C
Definition: External.C:196
Bool_t fUseBuiltinEventSelection
use builtin (old) event selection
TClonesArray * fTriggerPatchInfo
!trigger patch info array
Double_t fMaxPtHard
select maximum pt hard (MC)
Old centrality estimation (AliCentrality, works only on Run-1 PbPb and pPb)
TList * OpenFile(const char *fname)
Definition: DrawAnaELoss.C:65
Bool_t fNeedEmcalGeom
whether or not the task needs the emcal geometry
static Double_t GetParallelFraction(AliVParticle *part1, AliVParticle *part2)
EDataType_t fDataType
!data type (ESD or AOD)