AliPhysics  9b6b435 (9b6b435)
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 
64 
70  fForceBeamType(kNA),
71  fGeneralHistograms(kFALSE),
72  fCreateHisto(kTRUE),
73  fNeedEmcalGeom(kTRUE),
74  fCentBins(),
75  fCentralityEstimation(kNewCentrality),
76  fIsPythia(kFALSE),
77  fIsMonteCarlo(kFALSE),
78  fMCEventHeaderName(),
79  fCaloCellsName(),
80  fCaloTriggersName(),
81  fCaloTriggerPatchInfoName(),
82  fCentEst("V0M"),
83  fParticleCollArray(),
84  fClusterCollArray(),
85  fTriggerSelectionBitMap(0),
86  fMinCent(-1),
87  fMaxCent(-1),
88  fMinVz(-999),
89  fMaxVz(999),
90  fMaxVzDiff(-1),
91  fMinNVertCont(0),
92  fMinPtHard(-1),
93  fMaxPtHard(-1),
94  fMaxMinimumBiasPtHard(-1),
95  fAcceptedTriggerClasses(),
96  fRejectedTriggerClasses(),
97  fMCRejectFilter(kFALSE),
98  fPtHardAndJetPtFactor(0.),
99  fPtHardAndClusterPtFactor(0.),
100  fPtHardAndTrackPtFactor(0.),
101  fSwitchOffLHC15oFaultyBranches(kFALSE),
102  fEventSelectionAfterRun(kFALSE),
103  fUseAliEmcalList(kFALSE),
104  fUsePtHardBinScaling(kFALSE),
105  fSelectGeneratorName(),
106  fMinimumEventWeight(1e-6),
107  fMaximumEventWeight(1e6),
108  fInhibit(kFALSE),
109  fLocalInitialized(kFALSE),
110  fDataType(kAOD),
111  fGeom(0),
112  fCaloCells(0),
113  fCaloTriggers(0),
114  fTriggerPatchInfo(0),
115  fCent(-1),
116  fCentBin(-1),
117  fEPV0(-1.0),
118  fEPV0A(-1.0),
119  fEPV0C(-1.0),
120  fNVertCont(0),
121  fNVertSPDCont(0),
122  fFiredTriggerBitMap(0),
123  fFiredTriggerClasses(),
124  fBeamType(kNA),
125  fMCHeader(0),
126  fPythiaHeader(0),
127  fUseXsecFromHeader(false),
128  fPtHardBin(0),
129  fPtHard(0),
130  fNTrials(0),
131  fXsection(0),
132  fEventWeight(1),
133  fGeneratorName(),
134  fOutput(0),
135  fHistograms()
136 {
137  fVertex[0] = 0;
138  fVertex[1] = 0;
139  fVertex[2] = 0;
140  fVertexSPD[0] = 0;
141  fVertexSPD[1] = 0;
142  fVertexSPD[2] = 0;
143 }
144 
156  AliAnalysisTaskSE(name),
157  fForceBeamType(kNA),
158  fGeneralHistograms(kFALSE),
159  fCreateHisto(kTRUE),
160  fNeedEmcalGeom(kTRUE),
161  fCentBins(6),
162  fCentralityEstimation(kNewCentrality),
163  fIsPythia(kFALSE),
164  fIsMonteCarlo(kFALSE),
165  fMCEventHeaderName(),
166  fCaloCellsName(),
167  fCaloTriggersName(),
168  fCaloTriggerPatchInfoName(),
169  fCentEst("V0M"),
170  fParticleCollArray(),
171  fClusterCollArray(),
172  fTriggerSelectionBitMap(0),
173  fMinCent(-1),
174  fMaxCent(-1),
175  fMinVz(-999),
176  fMaxVz(999),
177  fMaxVzDiff(-1),
178  fMinNVertCont(0),
179  fMinPtHard(-1),
180  fMaxPtHard(-1),
181  fMaxMinimumBiasPtHard(-1),
182  fAcceptedTriggerClasses(),
183  fRejectedTriggerClasses(),
184  fMCRejectFilter(kFALSE),
185  fPtHardAndJetPtFactor(0.),
186  fPtHardAndClusterPtFactor(0.),
187  fPtHardAndTrackPtFactor(0.),
188  fSwitchOffLHC15oFaultyBranches(kFALSE),
189  fEventSelectionAfterRun(kFALSE),
190  fUseAliEmcalList(kFALSE),
191  fUsePtHardBinScaling(kFALSE),
192  fSelectGeneratorName(),
193  fMinimumEventWeight(1e-6),
194  fMaximumEventWeight(1e6),
195  fInhibit(kFALSE),
196  fLocalInitialized(kFALSE),
197  fDataType(kAOD),
198  fGeom(0),
199  fCaloCells(0),
200  fCaloTriggers(0),
201  fTriggerPatchInfo(0),
202  fCent(0),
203  fCentBin(-1),
204  fEPV0(-1.0),
205  fEPV0A(-1.0),
206  fEPV0C(-1.0),
207  fNVertCont(0),
208  fNVertSPDCont(0),
209  fFiredTriggerBitMap(0),
210  fFiredTriggerClasses(),
211  fBeamType(kNA),
212  fMCHeader(0),
213  fPythiaHeader(0),
214  fUseXsecFromHeader(false),
215  fPtHardBin(0),
216  fPtHard(0),
217  fNTrials(0),
218  fXsection(0),
219  fEventWeight(1),
220  fGeneratorName(),
221  fOutput(0),
222  fHistograms()
223 {
224  fVertex[0] = 0;
225  fVertex[1] = 0;
226  fVertex[2] = 0;
227  fVertexSPD[0] = 0;
228  fVertexSPD[1] = 0;
229  fVertexSPD[2] = 0;
230 
231  fCentBins[0] = 0;
232  fCentBins[1] = 10;
233  fCentBins[2] = 30;
234  fCentBins[3] = 50;
235  fCentBins[4] = 90;
236  fCentBins[5] = 100;
237 
238  if (fCreateHisto) DefineOutput(1, TList::Class());
239 }
240 
245 {
246  for (auto cont_it : fParticleCollArray) delete cont_it.second;
247  for (auto cont_it : fClusterCollArray) delete cont_it.second;
248 }
249 
270 {
271  if (fInhibit) {
272  AliWarningStream() << "The execution of this task is inhibited. Returning." << std::endl;
273  return;
274  }
275  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
276  if (mgr) {
277  AliVEventHandler *evhand = mgr->GetInputEventHandler();
278  if (evhand) {
279  if (evhand->InheritsFrom("AliESDInputHandler")) {
280  fDataType = kESD;
281  }
282  else {
283  fDataType = kAOD;
284  }
285  }
286  else {
287  AliError("Event handler not found!");
288  }
289  }
290  else {
291  AliError("Analysis manager not found!");
292  }
293 
294  if (!fCreateHisto)
295  return;
296 
297  OpenFile(1);
298  if(fUseAliEmcalList) {
299  auto emclist = new AliEmcalList;
300  if(fUsePtHardBinScaling) emclist->SetUseScaling(true);
301  emclist->SetNameXsec("fHistXsectionExternalFile");
302  emclist->SetNameTrials("fHistTrialsExternalFile");
303  fOutput = emclist;
304  } else {
305  fOutput = new TList();
306  }
307  fOutput->SetOwner(); // @suppress("Ambiguous problem")
308 
310 
311  if (!fGeneralHistograms) return;
312 
313  TH1* h = nullptr;
314 
315  if (fIsMonteCarlo) {
316  auto weight_bins = GenerateLogFixedBinArray(1000, fMinimumEventWeight, fMaximumEventWeight, true);
317 
318  h = new TH1F("fHistEventsVsPtHard", "fHistEventsVsPtHard", 1000, 0, 1000);
319  h->GetXaxis()->SetTitle("#it{p}_{T,hard} (GeV/#it{c})");
320  h->GetYaxis()->SetTitle("events");
321  fOutput->Add(h);
322  fHistograms["fHistEventsVsPtHard"] = h;
323 
324  h = new TH1F("fHistTrialsVsPtHard", "fHistTrialsVsPtHard", 1000, 0, 1000);
325  h->GetXaxis()->SetTitle("#it{p}_{T,hard} (GeV/#it{c})");
326  h->GetYaxis()->SetTitle("trials");
327  fOutput->Add(h);
328  fHistograms["fHistTrialsVsPtHard"] = h;
329 
330  h = new TProfile("fHistXsection", "fHistXsection", 50, 0, 50);
331  h->GetXaxis()->SetTitle("#it{p}_{T,hard} bin");
332  h->GetYaxis()->SetTitle("total integrated cross section (mb)");
333  fOutput->Add(h);
334  fHistograms["fHistXsection"] = h;
335 
336  h = new TH1F("fHistXsectionDistribution", "fHistXsectionDistribution", 1000, &weight_bins[0]);
337  h->GetXaxis()->SetTitle("total integrated cross section (mb)");
338  h->GetYaxis()->SetTitle("events");
339  fOutput->Add(h);
340  fHistograms["fHistXsectionDistribution"] = h;
341 
342  h = new TH1F("fHistEventWeights", "fHistEventWeights", 1000, &weight_bins[0]);
343  h->GetXaxis()->SetTitle("weight");
344  h->GetYaxis()->SetTitle("events");
345  fOutput->Add(h);
346  fHistograms["fHistEventWeights"] = h;
347 
348  h = new TH2F("fHistEventWeightsVsPtHard", "fHistEventWeightsVsPtHard", 1000, 0, 1000, 1000, &weight_bins[0]);
349  h->GetXaxis()->SetTitle("#it{p}_{T,hard} (GeV/#it{c})");
350  h->GetYaxis()->SetTitle("event weight");
351  fOutput->Add(h);
352  fHistograms["fHistEventWeightsVsPtHard"] = h;
353 
354  h = new TH1F("fHistEventsVsPtHardNoSel", "fHistEventsVsPtHardNoSel", 1000, 0, 1000);
355  h->GetXaxis()->SetTitle("#it{p}_{T,hard} (GeV/#it{c})");
356  h->GetYaxis()->SetTitle("events");
357  fOutput->Add(h);
358  fHistograms["fHistEventsVsPtHardNoSel"] = h;
359 
360  h = new TH1F("fHistTrialsVsPtHardNoSel", "fHistTrialsVsPtHardNoSel", 1000, 0, 1000);
361  h->GetXaxis()->SetTitle("#it{p}_{T,hard} (GeV/#it{c})");
362  h->GetYaxis()->SetTitle("trials");
363  fOutput->Add(h);
364  fHistograms["fHistTrialsVsPtHardNoSel"] = h;
365 
366  h = new TProfile("fHistXsectionNoSel", "fHistXsectionNoSel", 50, 0, 50);
367  h->GetXaxis()->SetTitle("#it{p}_{T,hard} bin");
368  h->GetYaxis()->SetTitle("total integrated cross section (mb)");
369  fOutput->Add(h);
370  fHistograms["fHistXsectionNoSel"] = h;
371 
372  h = new TH1F("fHistXsectionDistributionNoSel", "fHistXsectionDistributionNoSel", 1000, &weight_bins[0]);
373  h->GetXaxis()->SetTitle("total integrated cross section (mb)");
374  h->GetYaxis()->SetTitle("events");
375  fOutput->Add(h);
376  fHistograms["fHistXsectionDistributionNoSel"] = h;
377 
378  h = new TH1F("fHistEventWeightsNoSel", "fHistEventWeightsNoSel", 1000, &weight_bins[0]);
379  h->GetXaxis()->SetTitle("weight");
380  h->GetYaxis()->SetTitle("events");
381  fOutput->Add(h);
382  fHistograms["fHistEventWeightsNoSel"] = h;
383 
384  h = new TH2F("fHistEventWeightsVsPtHardNoSel", "fHistEventWeightsVsPtHardNoSel", 1000, 0, 1000, 1000, &weight_bins[0]);
385  h->GetXaxis()->SetTitle("#it{p}_{T,hard} (GeV/#it{c})");
386  h->GetYaxis()->SetTitle("event weight");
387  fOutput->Add(h);
388  fHistograms["fHistEventWeightsVsPtHardNoSel"] = h;
389 
390  h = new TH1F("fHistTrialsExternalFile", "fHistTrialsExternalFile", 50, 0, 50);
391  h->GetXaxis()->SetTitle("#it{p}_{T,hard} bin");
392  h->GetYaxis()->SetTitle("trials");
393  fOutput->Add(h);
394  fHistograms["fHistTrialsExternalFile"] = h;
395 
396  h = new TH1F("fHistEventsExternalFile", "fHistEventsExternalFile", 50, 0, 50);
397  h->GetXaxis()->SetTitle("#it{p}_{T,hard} bin");
398  h->GetYaxis()->SetTitle("total events");
399  fOutput->Add(h);
400  fHistograms["fHistEventsExternalFile"] = h;
401 
402  h = new TProfile("fHistXsectionExternalFile", "fHistXsectionExternalFile", 50, 0, 50);
403  h->GetXaxis()->SetTitle("#it{p}_{T,hard} bin");
404  h->GetYaxis()->SetTitle("total integrated cross section (mb)");
405  fOutput->Add(h);
406  fHistograms["fHistXsectionExternalFile"] = h;
407  }
408 
409  h = new TH1F("fHistZVertex","Z vertex position", 60, -30, 30);
410  h->GetXaxis()->SetTitle("V_{#it{z}}");
411  h->GetYaxis()->SetTitle("counts");
412  fOutput->Add(h);
413  fHistograms["fHistZVertex"] = h;
414 
415  h = new TH1F("fHistZVertexNoSel","Z vertex position (no event selection)", 60, -30, 30);
416  h->GetXaxis()->SetTitle("V_{#it{z}}");
417  h->GetYaxis()->SetTitle("counts");
418  fOutput->Add(h);
419  fHistograms["fHistZVertexNoSel"] = h;
420 
422  h = new TH1F("fHistCentrality","Event centrality distribution", 100, 0, 100);
423  h->GetXaxis()->SetTitle("Centrality (%)");
424  h->GetYaxis()->SetTitle("counts");
425  fOutput->Add(h);
426  fHistograms["fHistCentrality"] = h;
427 
428  h = new TH1F("fHistCentralityNoSel","Event centrality distribution (no event selection)", 100, 0, 100);
429  h->GetXaxis()->SetTitle("Centrality (%)");
430  h->GetYaxis()->SetTitle("counts");
431  fOutput->Add(h);
432  fHistograms["fHistCentralityNoSel"] = h;
433  }
434 
435  if (fForceBeamType != kpp) {
436  h = new TH1F("fHistEventPlane","Event plane", 120, -TMath::Pi(), TMath::Pi());
437  h->GetXaxis()->SetTitle("event plane");
438  h->GetYaxis()->SetTitle("counts");
439  fOutput->Add(h);
440  fHistograms["fHistEventPlane"] = h;
441 
442  h = new TH1F("fHistEventPlaneNoSel","Event plane (no event selection)", 120, -TMath::Pi(), TMath::Pi());
443  h->GetXaxis()->SetTitle("event plane");
444  h->GetYaxis()->SetTitle("counts");
445  fOutput->Add(h);
446  fHistograms["fHistEventPlaneNoSel"] = h;
447  }
448 
449  h = new TH1F("fHistEventRejection","Reasons to reject event",30,0,30);
450 #if ROOT_VERSION_CODE < ROOT_VERSION(6,4,2)
451  h->SetBit(TH1::kCanRebin);
452 #else
453  h->SetCanExtend(TH1::kAllAxes);
454 #endif
455  std::array<std::string, 10> labels = {"PhysSel", "Evt Gen Name", "Trg class (acc)", "Trg class (rej)", "Cent", "vertex contr.", "Vz", "VzSPD", "SelPtHardBin", "MCOutlier"};
456  int i = 1;
457  for (auto label : labels) {
458  h->GetXaxis()->SetBinLabel(i, label.c_str());
459  i++;
460  }
461  h->GetYaxis()->SetTitle("counts");
462  fOutput->Add(h);
463  fHistograms["fHistEventRejection"] = h;
464 
465  h = new TH1F("fHistTriggerClasses","fHistTriggerClasses",3,0,3);
466 #if ROOT_VERSION_CODE < ROOT_VERSION(6,4,2)
467  h->SetBit(TH1::kCanRebin);
468 #else
469  h->SetCanExtend(TH1::kAllAxes);
470 #endif
471  fOutput->Add(h);
472  fHistograms["fHistTriggerClasses"] = h;
473 
474  h = new TH1F("fHistTriggerClassesNoSel","fHistTriggerClassesNoSel",3,0,3);
475 #if ROOT_VERSION_CODE < ROOT_VERSION(6,4,2)
476  h->SetBit(TH1::kCanRebin);
477 #else
478  h->SetCanExtend(TH1::kAllAxes);
479 #endif
480  fOutput->Add(h);
481  fHistograms["fHistTriggerClassesNoSel"] = h;
482 
483  h = new TH1F("fHistEventCount","fHistEventCount",2,0,2);
484  h->GetXaxis()->SetBinLabel(1,"Accepted");
485  h->GetXaxis()->SetBinLabel(2,"Rejected");
486  h->GetYaxis()->SetTitle("counts");
487  fOutput->Add(h);
488  fHistograms["fHistEventCount"] = h;
489 
490  PostData(1, fOutput);
491 }
492 
508 {
509  if (eventSelected) {
510  if (fIsMonteCarlo) {
511  GetGeneralTH1("fHistEventsVsPtHard", true)->Fill(fPtHard);
512  GetGeneralTH1("fHistTrialsVsPtHard", true)->Fill(fPtHard, fNTrials);
513  GetGeneralTH1("fHistEventWeights", true)->Fill(fEventWeight);
514  GetGeneralTH2("fHistEventWeightsVsPtHard", true)->Fill(fPtHard, fEventWeight);
515  GetGeneralTH1("fHistXsectionDistribution", true)->Fill(fXsection);
516  GetGeneralTProfile("fHistXsection", true)->Fill(fPtHardBin, fXsection);
517  }
518 
519  GetGeneralTH1("fHistZVertex")->Fill(fVertex[2]);
520 
521  TH1* hCent = GetGeneralTH1("fHistCentrality");
522  if (hCent) hCent->Fill(fCent);
523 
524  TH1* hEventPlane = GetGeneralTH1("fHistEventPlane");
525  if (hEventPlane) hEventPlane->Fill(fEPV0);
526 
527  TH1* hTriggerClasses = GetGeneralTH1("fHistTriggerClasses");
528  for (auto fired_trg : fFiredTriggerClasses) hTriggerClasses->Fill(fired_trg.c_str(), 1);
529  }
530  else {
531  if (fIsMonteCarlo) {
532  GetGeneralTH1("fHistEventsVsPtHardNoSel", true)->Fill(fPtHard);
533  GetGeneralTH1("fHistTrialsVsPtHardNoSel", true)->Fill(fPtHard, fNTrials);
534  GetGeneralTH1("fHistEventWeightsNoSel", true)->Fill(fEventWeight);
535  GetGeneralTH2("fHistEventWeightsVsPtHardNoSel", true)->Fill(fPtHard, fEventWeight);
536  GetGeneralTH1("fHistXsectionDistributionNoSel", true)->Fill(fXsection);
537  GetGeneralTProfile("fHistXsectionNoSel", true)->Fill(fPtHardBin, fXsection);
538  }
539 
540  GetGeneralTH1("fHistZVertexNoSel", true)->Fill(fVertex[2]);
541 
542  TH1* hCent = GetGeneralTH1("fHistCentralityNoSel");
543  if (hCent) hCent->Fill(fCent);
544 
545  TH1* hEventPlane = GetGeneralTH1("fHistEventPlaneNoSel");
546  if (hEventPlane) hEventPlane->Fill(fEPV0);
547 
548  TH1* hTriggerClasses = GetGeneralTH1("fHistTriggerClassesNoSel", true);
549  for (auto fired_trg : fFiredTriggerClasses) hTriggerClasses->Fill(fired_trg.c_str(), 1);
550  }
551 
552  return kTRUE;
553 }
554 
575 {
576  if (fInhibit) {
577  AliWarningStream() << "The execution of this task is inhibited. Returning." << std::endl;
578  return;
579  }
580 
581  if (!fLocalInitialized) ExecOnce();
582 
583  if (!fLocalInitialized) return;
584 
585  if (!RetrieveEventObjects()) return;
586 
587  Bool_t eventSelected = IsEventSelected();
588 
590  if (eventSelected) {
591  GetGeneralTH1("fHistEventCount", true)->Fill("Accepted",1);
592  }
593  else {
594  GetGeneralTH1("fHistEventCount", true)->Fill("Rejected",1);
595  }
596 
597  FillGeneralHistograms(kFALSE);
598  if (eventSelected) FillGeneralHistograms(kTRUE);
599  }
600 
601  Bool_t runOk = kFALSE;
602  if (eventSelected || fEventSelectionAfterRun) runOk = Run();
603 
604  if (fCreateHisto && eventSelected && runOk) FillHistograms();
605 
606  if (fCreateHisto && fOutput) {
607  // information for this iteration of the UserExec in the container
608  PostData(1, fOutput);
609  }
610 }
611 
623 Bool_t AliAnalysisTaskEmcalLight::PythiaInfoFromFile(const char* currFile, Float_t &xsec, Float_t &trials, Int_t &pthard, Bool_t &useXsecFromHeader)
624 {
625 
626  TString file(currFile);
627  xsec = 0;
628  trials = 1;
629 
630  // Determine archive type
631  TString archivetype;
632  std::unique_ptr<TObjArray> walk(file.Tokenize("/"));
633  for(auto t : *walk){
634  TString &tok = static_cast<TObjString *>(t)->String();
635  if(tok.Contains(".zip")){
636  archivetype = tok;
637  Int_t pos = archivetype.Index(".zip");
638  archivetype.Replace(pos, archivetype.Length() - pos, "");
639  }
640  }
641  if(archivetype.Length()){
642  AliDebugStream(1) << "Auto-detected archive type " << archivetype << std::endl;
643  Ssiz_t pos1 = file.Index(archivetype,archivetype.Length(),0,TString::kExact);
644  Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
645  Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
646  file.Replace(pos+1,pos2-pos1,"");
647  } else {
648  // not an archive take the basename....
649  file.ReplaceAll(gSystem->BaseName(file.Data()),"");
650  }
651  AliDebugStream(1) << "File name: " << file << std::endl;
652 
653  // Build virtual file name
654  // Support for train tests
655  TString virtualFileName;
656  if(file.Contains("__alice")){
657  TString tmp(file);
658  Int_t pos = tmp.Index("__alice");
659  tmp.Replace(0, pos, "");
660  tmp.ReplaceAll("__", "/");
661  // cut out tag for archive and root file
662  // this needs a determin
663  std::unique_ptr<TObjArray> toks(tmp.Tokenize("/"));
664  TString tag = "_" + archivetype;
665  for(auto t : *toks){
666  TString &path = static_cast<TObjString *>(t)->String();
667  if(path.Contains(tag)){
668  Int_t posTag = path.Index(tag);
669  path.Replace(posTag, path.Length() - posTag, "");
670  }
671  virtualFileName += "/" + path;
672  }
673  } else {
674  virtualFileName = file;
675  }
676 
677  AliDebugStream(1) << "Physical file name " << file << ", virtual file name " << virtualFileName << std::endl;
678 
679  // Get the pt hard bin
680  TString strPthard(virtualFileName);
681 
682  /*
683  // Dead code - to be removed after testing phase
684  // Procedure will fail for everything else than the expected path name
685  strPthard.Remove(strPthard.Last('/'));
686  strPthard.Remove(strPthard.Last('/'));
687  if (strPthard.Contains("AOD")) strPthard.Remove(strPthard.Last('/'));
688  strPthard.Remove(0,strPthard.Last('/')+1);
689  if (strPthard.IsDec()) pthard = strPthard.Atoi();
690  else
691  AliWarningStream() << "Could not extract file number from path " << strPthard << std::endl;
692  */
693 
694  // New implementation : pattern matching
695  // Reason: Implementation valid only for old productions (new productions swap run number and pt-hard bin)
696  // Idea: Don't use the position in the string but the match different informations
697  // + Year clearly 2000+
698  // + Run number can be match to the one in the event
699  // + If we know it is not year or run number, it must be the pt-hard bin if we start from the beginning
700  // The procedure is only valid for the current implementations and unable to detect non-pt-hard bins
701  // It will also fail in case of arbitrary file names
702 
703  bool binfound = false;
704  std::unique_ptr<TObjArray> tokens(strPthard.Tokenize("/"));
705  for(auto t : *tokens) {
706  TString &tok = static_cast<TObjString *>(t)->String();
707  if(tok.IsDec()){
708  Int_t number = tok.Atoi();
709  if(number > 2000 && number < 3000){
710  // Year
711  continue;
712  } else if(number == fInputHandler->GetEvent()->GetRunNumber()){
713  // Run number
714  continue;
715  } else {
716  if(!binfound){
717  // the first number that is not one of the two must be the pt-hard bin
718  binfound = true;
719  pthard = number;
720  break;
721  }
722  }
723  }
724  }
725  if(!binfound) {
726  AliErrorStream() << "Could not extract file number from path " << strPthard << std::endl;
727  } else {
728  AliInfoStream() << "Auto-detecting pt-hard bin " << pthard << std::endl;
729  }
730 
731  AliInfoStream() << "File: " << file << std::endl;
732 
733  // 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
734  std::unique_ptr<TFile> fxsec(TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")));
735 
736  if (!fxsec) {
737  // next trial fetch the histgram file
738  fxsec = std::unique_ptr<TFile>(TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root")));
739  if (!fxsec){
740  AliErrorStream() << "Failed reading cross section from file " << file << std::endl;
741  useXsecFromHeader = true;
742  return kFALSE; // not a severe condition but inciate that we have no information
743  }
744  else {
745  // find the tlist we want to be independtent of the name so use the Tkey
746  TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
747  if (!key) return kFALSE;
748  TList *list = dynamic_cast<TList*>(key->ReadObj());
749  if (!list) return kFALSE;
750  TProfile *xSecHist = static_cast<TProfile*>(list->FindObject("h1Xsec"));
751  // check for failure
752  if(!xSecHist->GetEntries()) {
753  // No cross seciton information available - fall back to raw
754  AliErrorStream() << "No cross section information available in file " << fxsec->GetName() <<" - fall back to cross section in PYTHIA header" << std::endl;
755  useXsecFromHeader = true;
756  } else {
757  // Cross section histogram filled - take it from there
758  xsec = xSecHist->GetBinContent(1);
759  if(!xsec) AliErrorStream() << GetName() << ": Cross section 0 for file " << file << std::endl;
760  useXsecFromHeader = false;
761  }
762  trials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
763  }
764  } else { // no tree pyxsec.root
765  TTree *xtree = (TTree*)fxsec->Get("Xsection");
766  if (!xtree) return kFALSE;
767  UInt_t ntrials = 0;
768  Double_t xsection = 0;
769  xtree->SetBranchAddress("xsection",&xsection);
770  xtree->SetBranchAddress("ntrials",&ntrials);
771  xtree->GetEntry(0);
772  trials = ntrials;
773  xsec = xsection;
774  }
775  return kTRUE;
776 }
777 
792 {
793  if (!fIsPythia) return kTRUE;
794 
795  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
796  if (!tree) {
797  AliError(Form("%s - UserNotify: No current tree!",GetName()));
798  return kFALSE;
799  }
800 
801  Float_t xsection = 0;
802  Float_t trials = 0;
803  Int_t pthardbin = 0;
804  Bool_t useXsecFromHeader = false;
805 
806  TFile *curfile = tree->GetCurrentFile();
807  if (!curfile) {
808  AliError(Form("%s - UserNotify: No current file!",GetName()));
809  return kFALSE;
810  }
811 
812  TChain *chain = dynamic_cast<TChain*>(tree);
813  if (chain) tree = chain->GetTree();
814 
815  Int_t nevents = tree->GetEntriesFast();
816 
817  Bool_t res = PythiaInfoFromFile(curfile->GetName(), xsection, trials, pthardbin, useXsecFromHeader);
818 
819  fPtHardBin = pthardbin >= 0 ? pthardbin : 0;
820  fUseXsecFromHeader = useXsecFromHeader;
821 
822  if (!res) return kTRUE;
823 
825  GetGeneralTH1("fHistTrialsExternalFile", true)->Fill(fPtHardBin, trials);
826  if(!useXsecFromHeader) GetGeneralTProfile("fHistXsectionExternalFile", true)->Fill(fPtHardBin, xsection);
827  GetGeneralTH1("fHistEventsExternalFile", true)->Fill(fPtHardBin, nevents);
828  }
829 
830  return kTRUE;
831 }
832 
844 {
845  if (!InputEvent()) {
846  AliError(Form("%s: Could not retrieve event! Returning!", GetName()));
847  return;
848  }
849 
850  if (fNeedEmcalGeom && !fGeom) {
851  fGeom = AliEMCALGeometry::GetInstanceFromRunNumber(InputEvent()->GetRunNumber());
852  if (!fGeom) {
853  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()));
854  return;
855  }
856  }
857 
858  if (fSwitchOffLHC15oFaultyBranches && dynamic_cast<AliAODEvent*>(InputEvent())) {
859  TTree *aodTree = AliAnalysisManager::GetAnalysisManager()->GetTree();
860  aodTree->SetBranchStatus("D0toKpi.fPx", 0);
861  aodTree->SetBranchStatus("D0toKpi.fPy", 0);
862  aodTree->SetBranchStatus("D0toKpi.fPz", 0);
863  aodTree->SetBranchStatus("D0toKpi.fd0", 0);
864  aodTree->SetBranchStatus("Charm3Prong.fPx", 0);
865  aodTree->SetBranchStatus("Charm3Prong.fPy", 0);
866  aodTree->SetBranchStatus("Charm3Prong.fPz", 0);
867  aodTree->SetBranchStatus("Charm3Prong.fd0", 0);
868  aodTree->SetBranchStatus("Dstar.fPx", 0);
869  aodTree->SetBranchStatus("Dstar.fPy", 0);
870  aodTree->SetBranchStatus("Dstar.fPz", 0);
871  aodTree->SetBranchStatus("Dstar.fd0", 0);
872  }
873 
874  //Load all requested track branches - each container knows name already
875  for (auto cont_it : fParticleCollArray) cont_it.second->SetArray(InputEvent());
876 
877  //Load all requested cluster branches - each container knows name already
878  for (auto cont_it : fClusterCollArray) cont_it.second->SetArray(InputEvent());
879 
880  if (!fCaloCellsName.IsNull() && !fCaloCells) {
881  fCaloCells = dynamic_cast<AliVCaloCells*>(InputEvent()->FindListObject(fCaloCellsName));
882  if (!fCaloCells) {
883  AliError(Form("%s: Could not retrieve cells %s!", GetName(), fCaloCellsName.Data()));
884  return;
885  }
886  }
887 
888  if (!fCaloTriggersName.IsNull() && !fCaloTriggers) {
889  fCaloTriggers = dynamic_cast<AliVCaloTrigger*>(InputEvent()->FindListObject(fCaloTriggersName));
890  if (!fCaloTriggers) {
891  AliError(Form("%s: Could not retrieve calo triggers %s!", GetName(), fCaloTriggersName.Data()));
892  return;
893  }
894  }
895 
896  if (!fCaloTriggerPatchInfoName.IsNull() && !fTriggerPatchInfo) {
897  fTriggerPatchInfo = GetArrayFromEvent(fCaloTriggerPatchInfoName.Data(),"AliEMCALTriggerPatchInfo");
898  if (!fTriggerPatchInfo) {
899  AliError(Form("%s: Could not retrieve calo trigger patch info %s!", GetName(), fCaloTriggerPatchInfoName.Data()));
900  return;
901  }
902 
903  }
904 
905  fLocalInitialized = kTRUE;
906 }
907 
914 {
915  if (fForceBeamType != kNA)
916  return fForceBeamType;
917 
918  AliESDEvent *esd = dynamic_cast<AliESDEvent*>(InputEvent());
919  if (esd) {
920  const AliESDRun *run = esd->GetESDRun();
921  TString beamType = run->GetBeamType();
922  if (beamType == "p-p")
923  return kpp;
924  else if (beamType == "A-A")
925  return kAA;
926  else if (beamType == "p-A")
927  return kpA;
928  else
929  return kNA;
930  } else {
931  Int_t runNumber = InputEvent()->GetRunNumber();
932  // All run number ranges taken from the RCT
933  if ((runNumber >= 136833 && runNumber <= 139517) || // LHC10h
934  (runNumber >= 167693 && runNumber <= 170593) || // LHC11h
935  (runNumber >= 244824 && runNumber <= 246994)) { // LHC15o
936  return kAA;
937  } else if ((runNumber >= 188356 && runNumber <= 188366) || // LHC12g
938  (runNumber >= 195164 && runNumber <= 197388) || // LHC13b-f
939  (runNumber >= 265015 && runNumber <= 267166)) { // LHC16q-t
940  return kpA;
941  } else {
942  return kpp;
943  }
944  }
945 }
946 
969 {
970  TH1* hEventRejection = GetGeneralTH1("fHistEventRejection", true);
971 
973  if (fGeneralHistograms) hEventRejection->Fill("PhysSel",1);
974  return kFALSE;
975  }
976 
977  if (!fSelectGeneratorName.IsNull() && !fGeneratorName.IsNull()) {
978  if (!fGeneratorName.Contains(fSelectGeneratorName)) {
979  if (fGeneralHistograms) hEventRejection->Fill("Evt Gen Name",1);
980  return kFALSE;
981  }
982  }
983 
984  Bool_t acceptedTrgClassFound = kFALSE;
985  if (fAcceptedTriggerClasses.size() > 0) {
986  for (auto acc_trg : fAcceptedTriggerClasses) {
987  std::string teststring(acc_trg);
988  bool fullmatch(false);
989  auto posexact = acc_trg.find("EXACT");
990  if(posexact != std::string::npos) {
991  fullmatch = true;
992  teststring.erase(posexact, 5);
993  }
994  for (auto fired_trg : fFiredTriggerClasses) {
995  bool classmatch = fullmatch ? teststring == fired_trg : fired_trg.find(teststring) != std::string::npos;
996  if (classmatch) {
997  acceptedTrgClassFound = kTRUE;
998  break;
999  }
1000  }
1001  if (acceptedTrgClassFound) break;
1002  }
1003 
1004  if (!acceptedTrgClassFound) {
1005  if (fGeneralHistograms) hEventRejection->Fill("Trg class (acc)",1);
1006  return kFALSE;
1007  }
1008  }
1009 
1010  if (fRejectedTriggerClasses.size() > 0) {
1011  for (auto rej_trg : fRejectedTriggerClasses) {
1012  std::string teststring(rej_trg);
1013  bool fullmatch(false);
1014  auto posexact = rej_trg.find("EXACT");
1015  if(posexact != std::string::npos) {
1016  fullmatch = true;
1017  teststring.erase(posexact, 5);
1018  }
1019  for (auto fired_trg : fFiredTriggerClasses) {
1020  bool classmatch = fullmatch ? teststring == fired_trg : fired_trg.find(teststring) != std::string::npos;
1021  if (classmatch) {
1022  if (fGeneralHistograms) hEventRejection->Fill("Trg class (rej)",1);
1023  return kFALSE;
1024  }
1025  }
1026  }
1027  }
1028 
1029  if (fMinCent < fMaxCent && fMaxCent > 0) {
1030  if (fCent < fMinCent || fCent > fMaxCent) {
1031  if (fGeneralHistograms) hEventRejection->Fill("Cent",1);
1032  return kFALSE;
1033  }
1034  }
1035 
1036  if (fNVertCont < fMinNVertCont) {
1037  if (fGeneralHistograms) hEventRejection->Fill("vertex contr.",1);
1038  return kFALSE;
1039  }
1040 
1041  if (fMinVz < fMaxVz) {
1042  if (fVertex[2] < fMinVz || fVertex[2] > fMaxVz) {
1043  if (fGeneralHistograms) hEventRejection->Fill("Vz",1);
1044  return kFALSE;
1045  }
1046  }
1047 
1048  if (fMaxVzDiff >= 0) {
1049  if (fNVertSPDCont > 0) {
1050  Double_t vzSPD = fVertexSPD[2];
1051  Double_t dvertex = TMath::Abs(fVertex[2] - vzSPD);
1052  //if difference larger than fZvertexDiff
1053  if (dvertex > fMaxVzDiff) {
1054  if (fGeneralHistograms) hEventRejection->Fill("VzSPD",1);
1055  return kFALSE;
1056  }
1057  }
1058  }
1059 
1060  if (fMinPtHard >= 0 && fPtHard < fMinPtHard) {
1061  if (fGeneralHistograms) hEventRejection->Fill("SelPtHardBin",1);
1062  return kFALSE;
1063  }
1064 
1065  if (fMaxPtHard >= 0 && fPtHard >= fMaxPtHard) {
1066  if (fGeneralHistograms) hEventRejection->Fill("SelPtHardBin",1);
1067  return kFALSE;
1068  }
1069 
1071  if (fGeneralHistograms) hEventRejection->Fill("SelPtHardBin",1);
1072  return kFALSE;
1073  }
1074 
1075  // Reject filter for MC data
1076  if (!CheckMCOutliers()) {
1077  if (fGeneralHistograms) hEventRejection->Fill("MCOutlier",1);
1078  return kFALSE;
1079  }
1080 
1081  return kTRUE;
1082 }
1083 
1092 TClonesArray *AliAnalysisTaskEmcalLight::GetArrayFromEvent(const char *name, const char *clname)
1093 {
1094  TClonesArray *arr = 0;
1095  TString sname(name);
1096  if (!sname.IsNull()) {
1097  arr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(sname));
1098  if (!arr) {
1099  AliWarning(Form("%s: Could not retrieve array with name %s!", GetName(), name));
1100  return 0;
1101  }
1102  } else {
1103  return 0;
1104  }
1105 
1106  if (!clname)
1107  return arr;
1108 
1109  TString objname(arr->GetClass()->GetName());
1110  TClass cls(objname);
1111  if (!cls.InheritsFrom(clname)) {
1112  AliWarning(Form("%s: Objects of type %s in %s are not inherited from %s!",
1113  GetName(), cls.GetName(), name, clname));
1114  return 0;
1115  }
1116  return arr;
1117 }
1118 
1124 {
1125  fVertex[0] = 0;
1126  fVertex[1] = 0;
1127  fVertex[2] = 0;
1128  fNVertCont = 0;
1129 
1130  fVertexSPD[0] = 0;
1131  fVertexSPD[1] = 0;
1132  fVertexSPD[2] = 0;
1133  fNVertSPDCont = 0;
1134 
1135  fFiredTriggerClasses.clear();
1136  std::stringstream firedClasses(InputEvent()->GetFiredTriggerClasses().Data());
1137  while (firedClasses.good()) {
1138  std::string trgClass;
1139  firedClasses >> trgClass;
1140  if (!trgClass.empty()) fFiredTriggerClasses.push_back(trgClass);
1141  }
1142 
1143  if (fDataType == kESD) {
1144  fFiredTriggerBitMap = static_cast<AliInputEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler())->IsEventSelected();
1145  }
1146  else {
1147  fFiredTriggerBitMap = static_cast<AliVAODHeader*>(InputEvent()->GetHeader())->GetOfflineTrigger();
1148  }
1149 
1150  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
1151  if (vert) {
1152  vert->GetXYZ(fVertex);
1153  fNVertCont = vert->GetNContributors();
1154  }
1155 
1156  const AliVVertex *vertSPD = InputEvent()->GetPrimaryVertexSPD();
1157  if (vertSPD) {
1158  vertSPD->GetXYZ(fVertexSPD);
1159  fNVertSPDCont = vertSPD->GetNContributors();
1160  }
1161 
1162  fBeamType = GetBeamType();
1163 
1164  fCent = 99;
1165  fCentBin = -1;
1166  fEPV0 = -999;
1167  fEPV0A = -999;
1168  fEPV0C = -999;
1169 
1171  // New centrality estimation (AliMultSelection)
1172  // See https://twiki.cern.ch/twiki/bin/viewauth/ALICE/AliMultSelectionCalibStatus for calibration status period-by-period)
1173  AliMultSelection *MultSelection = static_cast<AliMultSelection*>(InputEvent()->FindListObject("MultSelection"));
1174  if (MultSelection) {
1175  fCent = MultSelection->GetMultiplicityPercentile(fCentEst.Data());
1176  }
1177  else {
1178  AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
1179  }
1180  }
1181  else if (fCentralityEstimation == kOldCentrality) {
1182  // Old centrality estimation (AliCentrality, works only on Run-1 PbPb and pPb)
1183  AliCentrality *aliCent = InputEvent()->GetCentrality();
1184  if (aliCent) {
1185  fCent = aliCent->GetCentralityPercentile(fCentEst.Data());
1186  }
1187  else {
1188  AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
1189  }
1190  }
1191  if (!fCentBins.empty() && fCentralityEstimation != kNoCentrality) {
1192  for (auto cent_it = fCentBins.begin(); cent_it != fCentBins.end() - 1; cent_it++) {
1193  if (fCent >= *cent_it && fCent < *(cent_it+1)) fCentBin = cent_it - fCentBins.begin();
1194  }
1195  }
1196  else {
1197  fCentBin = 0;
1198  }
1199 
1200  if (fBeamType == kAA || fBeamType == kpA ) {
1201  AliEventplane *aliEP = InputEvent()->GetEventplane();
1202  if (aliEP) {
1203  fEPV0 = aliEP->GetEventplane("V0" ,InputEvent());
1204  fEPV0A = aliEP->GetEventplane("V0A",InputEvent());
1205  fEPV0C = aliEP->GetEventplane("V0C",InputEvent());
1206  } else {
1207  AliWarning(Form("%s: Could not retrieve event plane information!", GetName()));
1208  }
1209  }
1210 
1211  if (fIsMonteCarlo && MCEvent()) {
1212  AliGenEventHeader* header = MCEvent()->GenEventHeader();
1213  if (fMCEventHeaderName.IsNull()) {
1214  fMCHeader = header;
1215  }
1216  else {
1217  if (header->InheritsFrom(fMCEventHeaderName)) {
1218  fMCHeader = header;
1219  }
1220  else if (header->InheritsFrom("AliGenCocktailEventHeader")) {
1221  AliGenCocktailEventHeader* cocktailHeader = static_cast<AliGenCocktailEventHeader*>(header);
1222  TList* headers = cocktailHeader->GetHeaders();
1223  for (auto obj : *headers) { // @suppress("Symbol is not resolved")
1224  if (obj->InheritsFrom(fMCEventHeaderName)){
1225  fMCHeader = static_cast<AliGenEventHeader*>(obj);
1226  break;
1227  }
1228  }
1229  }
1230  }
1231  if (fMCHeader) {
1232  fEventWeight = fMCHeader->EventWeight();
1233  if (fIsPythia) {
1234  fPythiaHeader = static_cast<AliGenPythiaEventHeader*>(fMCHeader);
1235  fPtHard = fPythiaHeader->GetPtHard();
1236  fXsection = fPythiaHeader->GetXsection();
1237  fNTrials = fPythiaHeader->Trials();
1238  if(fUseXsecFromHeader) GetGeneralTProfile("fHistXsectionExternalFile", true)->Fill(fPtHardBin, fXsection);
1239  }
1240  }
1241  }
1242 
1243  for (auto cont_it : fParticleCollArray) cont_it.second->NextEvent(InputEvent());
1244  for (auto cont_it : fClusterCollArray) cont_it.second->NextEvent(InputEvent());
1245 
1246  return kTRUE;
1247 }
1248 
1257 AliParticleContainer* AliAnalysisTaskEmcalLight::AddParticleContainer(std::string branchName, std::string contName)
1258 {
1259  if (branchName.size() == 0) return 0;
1260 
1261  AliParticleContainer* cont = 0;
1262 
1263  if (branchName == "tracks" || branchName == "Tracks") cont = new AliTrackContainer(branchName.c_str());
1264  else if (branchName == "mcparticles") cont = new AliMCParticleContainer(branchName.c_str());
1265  else cont = new AliParticleContainer(branchName.c_str());
1266 
1267  if (contName.size() > 0) cont->SetName(contName.c_str());
1268 
1269  AdoptParticleContainer(cont);
1270 
1271  return cont;
1272 }
1273 
1282 AliClusterContainer* AliAnalysisTaskEmcalLight::AddClusterContainer(std::string branchName, std::string contName)
1283 {
1284  if (branchName.size() == 0) return 0;
1285 
1286  AliClusterContainer* cont = new AliClusterContainer(branchName.c_str());
1287 
1288  if (contName.size() > 0) cont->SetName(contName.c_str());
1289 
1290  AdoptClusterContainer(cont);
1291 
1292  return cont;
1293 }
1294 
1301 {
1302  std::map<std::string, AliParticleContainer*>::const_iterator cont_it = fParticleCollArray.find(name);
1303  if (cont_it != fParticleCollArray.end()) return cont_it->second;
1304  else return nullptr;
1305 }
1306 
1313 {
1314  std::map<std::string, AliClusterContainer*>::const_iterator cont_it = fClusterCollArray.find(name);
1315  if (cont_it != fClusterCollArray.end()) return cont_it->second;
1316  else return nullptr;
1317 }
1318 
1325 {
1326  if (!(InputEvent()->FindListObject(obj->GetName()))) {
1327  InputEvent()->AddObject(obj);
1328  }
1329  else {
1330  if (!attempt) {
1331  AliFatal(Form("%s: Container with name %s already present. Aborting", GetName(), obj->GetName()));
1332  }
1333  }
1334 }
1335 
1344 {
1345 
1346  if (!fGeom) {
1347  AliWarning(Form("%s - AliAnalysisTaskEmcalBase::IsTrackInEmcalAcceptance - Geometry is not available!", GetName()));
1348  return kFALSE;
1349  }
1350 
1351  Double_t minPhi = fGeom->GetArm1PhiMin() - edges;
1352  Double_t maxPhi = fGeom->GetArm1PhiMax() + edges;
1353 
1354  if (part->Phi() > minPhi && part->Phi() < maxPhi) {
1355  return kTRUE;
1356  }
1357  else {
1358  return kFALSE;
1359  }
1360 }
1361 
1363 {
1364  axis->SetBinLabel(1, "NullObject");
1365  axis->SetBinLabel(2, "Pt");
1366  axis->SetBinLabel(3, "Acceptance");
1367  axis->SetBinLabel(4, "MCLabel");
1368  axis->SetBinLabel(5, "BitMap");
1369  axis->SetBinLabel(6, "HF cut");
1370  axis->SetBinLabel(7, "Bit6");
1371  axis->SetBinLabel(8, "NotHybridTrack");
1372  axis->SetBinLabel(9, "MCFlag");
1373  axis->SetBinLabel(10, "MCGenerator");
1374  axis->SetBinLabel(11, "ChargeCut");
1375  axis->SetBinLabel(12, "MinDistanceTPCSectorEdge");
1376  axis->SetBinLabel(13, "Bit12");
1377  axis->SetBinLabel(14, "IsEMCal");
1378  axis->SetBinLabel(15, "Time");
1379  axis->SetBinLabel(16, "Energy");
1380  axis->SetBinLabel(17, "ExoticCut");
1381  axis->SetBinLabel(18, "Bit17");
1382  axis->SetBinLabel(19, "Area");
1383  axis->SetBinLabel(20, "AreaEmc");
1384  axis->SetBinLabel(21, "ZLeadingCh");
1385  axis->SetBinLabel(22, "ZLeadingEmc");
1386  axis->SetBinLabel(23, "NEF");
1387  axis->SetBinLabel(24, "MinLeadPt");
1388  axis->SetBinLabel(25, "MaxTrackPt");
1389  axis->SetBinLabel(26, "MaxClusterPt");
1390  axis->SetBinLabel(27, "Flavour");
1391  axis->SetBinLabel(28, "TagStatus");
1392  axis->SetBinLabel(29, "MinNConstituents");
1393  axis->SetBinLabel(30, "Bit29");
1394  axis->SetBinLabel(31, "Bit30");
1395  axis->SetBinLabel(32, "Bit31");
1396 }
1397 
1404 Double_t AliAnalysisTaskEmcalLight::GetParallelFraction(AliVParticle* part1, AliVParticle* part2)
1405 {
1406  TVector3 vect1(part1->Px(), part1->Py(), part1->Pz());
1407  TVector3 vect2(part2->Px(), part2->Py(), part2->Pz());
1408  Double_t z = (vect1 * vect2) / (vect2 * vect2);
1409  return z;
1410 }
1411 
1418 Double_t AliAnalysisTaskEmcalLight::GetParallelFraction(const TVector3& vect1, AliVParticle* part2)
1419 {
1420  TVector3 vect2(part2->Px(), part2->Py(), part2->Pz());
1421  Double_t z = (vect1 * vect2) / (vect2 * vect2);
1422  return z;
1423 }
1424 
1433 void AliAnalysisTaskEmcalLight::GetEtaPhiDiff(const AliVTrack *t, const AliVCluster *v, Double_t &phidiff, Double_t &etadiff)
1434 {
1435  phidiff = 999;
1436  etadiff = 999;
1437 
1438  if (!t||!v) return;
1439 
1440  Double_t veta = t->GetTrackEtaOnEMCal();
1441  Double_t vphi = t->GetTrackPhiOnEMCal();
1442 
1443  Float_t pos[3] = {0};
1444  v->GetPosition(pos);
1445  TVector3 cpos(pos);
1446  Double_t ceta = cpos.Eta();
1447  Double_t cphi = cpos.Phi();
1448  etadiff=veta-ceta;
1449  phidiff=TVector2::Phi_mpi_pi(vphi-cphi);
1450 }
1451 
1458 {
1459  Byte_t ret = 0;
1460  if (t->TestBit(BIT(22)) && !t->TestBit(BIT(23)))
1461  ret = 1;
1462  else if (!t->TestBit(BIT(22)) && t->TestBit(BIT(23)))
1463  ret = 2;
1464  else if (t->TestBit(BIT(22)) && t->TestBit(BIT(23)))
1465  ret = 3;
1466  return ret;
1467 }
1468 
1478 Byte_t AliAnalysisTaskEmcalLight::GetTrackType(const AliAODTrack *aodTrack, UInt_t filterBit1, UInt_t filterBit2)
1479 {
1480 
1481  Int_t res = 0;
1482 
1483  if (aodTrack->TestFilterBit(filterBit1)) {
1484  res = 0;
1485  }
1486  else if (aodTrack->TestFilterBit(filterBit2)) {
1487  if ((aodTrack->GetStatus()&AliVTrack::kITSrefit)!=0) {
1488  res = 1;
1489  }
1490  else {
1491  res = 2;
1492  }
1493  }
1494  else {
1495  res = 3;
1496  }
1497 
1498  return res;
1499 }
1500 
1507 {
1508  EBeamType_t b = kpp;
1509  if ((runnumber >= 136833 && runnumber <= 139517) || // LHC10h Run-1 (Pb-Pb)
1510  (runnumber >= 167693 && runnumber <= 170593) || // LHC11h Run-1 (Pb-Pb)
1511  (runnumber >= 244824 && runnumber <= 246994) //|| // LHC15o Run-2 (Pb-Pb)
1512  //(runnumber >= 295581 && runnumber <= 297624) // LHC18q+r Run-2 (Pb-Pb)
1513  )
1514  {
1515  b = kAA;
1516  }
1517  else if ((runnumber > 188356 && runnumber <= 188503) || // LHC12g Run-1 (p-Pb pilot)
1518  (runnumber >= 195164 && runnumber <= 197388) || // LHC13b,c,d,e,f Run-1 (p-Pb)
1519  (runnumber >= 265077 && runnumber <= 267166)) { // LHC16 Run-2 (p-Pb)
1520  b = kpA;
1521  }
1522  return b;
1523 }
1524 
1531 {
1532  if (!fPythiaHeader || !fMCRejectFilter) return kTRUE;
1533 
1534  // Condition 1: Pythia jet / pT-hard > factor
1535  if (fPtHardAndJetPtFactor > 0.) {
1536  AliTLorentzVector jet;
1537 
1538  Int_t nTriggerJets = fPythiaHeader->NTriggerJets();
1539 
1540  AliDebug(1,Form("Njets: %d, pT Hard %f",nTriggerJets, fPtHard));
1541 
1542  Float_t tmpjet[]={0,0,0,0};
1543  for (Int_t ijet = 0; ijet< nTriggerJets; ijet++) {
1544  fPythiaHeader->TriggerJet(ijet, tmpjet);
1545 
1546  jet.SetPxPyPzE(tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3]);
1547 
1548  AliDebug(1,Form("jet %d; pycell jet pT %f",ijet, jet.Pt()));
1549 
1550  //Compare jet pT and pt Hard
1551  if (jet.Pt() > fPtHardAndJetPtFactor * fPtHard) {
1552  AliInfo(Form("Reject jet event with : pT Hard %2.2f, pycell jet pT %2.2f, rejection factor %1.1f\n", fPtHard, jet.Pt(), fPtHardAndJetPtFactor));
1553  return kFALSE;
1554  }
1555  }
1556  }
1557  // end condition 1
1558 
1559  // Condition 2 : Reconstructed EMCal cluster pT / pT-hard > factor
1560  if (fPtHardAndClusterPtFactor > 0.) {
1561  AliClusterContainer* mccluscont = GetClusterContainer(0);
1562  if ((Bool_t)mccluscont) {
1563  for (auto cluster : mccluscont->all()) {// Not cuts applied ; use accept for cuts
1564  Float_t ecluster = cluster->E();
1565 
1566  if (ecluster > (fPtHardAndClusterPtFactor * fPtHard)) {
1567  AliInfo(Form("Reject : ecluster %2.2f, calo %d, factor %2.2f, ptHard %f",ecluster,cluster->GetType(),fPtHardAndClusterPtFactor,fPtHard));
1568  return kFALSE;
1569  }
1570  }
1571  }
1572  }
1573  // end condition 2
1574 
1575  // condition 3 : Reconstructed track pT / pT-hard >factor
1576  if (fPtHardAndTrackPtFactor > 0.) {
1577  AliMCParticleContainer* mcpartcont = dynamic_cast<AliMCParticleContainer*>(GetParticleContainer(0));
1578  if ((Bool_t)mcpartcont) {
1579  for (auto mctrack : mcpartcont->all()) {// Not cuts applied ; use accept for cuts
1580  Float_t trackpt = mctrack->Pt();
1581  if (trackpt > (fPtHardAndTrackPtFactor * fPtHard) ) {
1582  AliInfo(Form("Reject : track %2.2f, factor %2.2f, ptHard %f", trackpt, fPtHardAndTrackPtFactor, fPtHard));
1583  return kFALSE;
1584  }
1585  }
1586  }
1587  }
1588  // end condition 3
1589 
1590  return kTRUE;
1591 }
1592 
1602 {
1603  Double_t dphi = -999;
1604  const Double_t tpi = TMath::TwoPi();
1605 
1606  if (phia < 0) phia += tpi;
1607  else if (phia > tpi) phia -= tpi;
1608  if (phib < 0) phib += tpi;
1609  else if (phib > tpi) phib -= tpi;
1610  dphi = phib - phia;
1611  if (dphi < rangeMin) dphi += tpi;
1612  else if (dphi > rangeMax) dphi -= tpi;
1613 
1614  return dphi;
1615 }
1616 
1626 void AliAnalysisTaskEmcalLight::GenerateFixedBinArray(int n, double min, double max, std::vector<double>& array, bool last)
1627 {
1628  double binWidth = (max - min) / n;
1629  double v = min;
1630  if (last) n++;
1631  for (int i = 0; i < n; i++) {
1632  array.push_back(v);
1633  v += binWidth;
1634  }
1635 }
1636 
1646 std::vector<double> AliAnalysisTaskEmcalLight::GenerateFixedBinArray(int n, double min, double max, bool last)
1647 {
1648  std::vector<double> array;
1649  GenerateFixedBinArray(n, min, max, array, last);
1650  return array;
1651 }
1652 
1662 void AliAnalysisTaskEmcalLight::GenerateLogFixedBinArray(int n, double min, double max, std::vector<double>& array, bool last)
1663 {
1664  if (min <= 0 || max < min) {
1665  AliErrorClassStream() << "Cannot generate a log scale fixed-bin array with limits " << min << ", " << max << std::endl;
1666  return;
1667  }
1668  double binWidth = std::pow(max / min, 1.0 / n);
1669  double v = min;
1670  if (last) n++;
1671  for (int i = 0; i < n; i++) {
1672  array.push_back(v);
1673  v *= binWidth;
1674  }
1675 }
1676 
1686 std::vector<double> AliAnalysisTaskEmcalLight::GenerateLogFixedBinArray(int n, double min, double max, bool last)
1687 {
1688  std::vector<double> array;
1689  GenerateLogFixedBinArray(n, min, max, array, last);
1690  return array;
1691 }
1692 
1693 
1694 TH1* AliAnalysisTaskEmcalLight::GetGeneralTH1(const char* name, bool warn)
1695 {
1696  auto search = fHistograms.find(name);
1697  if (search != fHistograms.end()) {
1698  return search->second;
1699  }
1700  else {
1701  if (warn) AliErrorStream() << "Could not find histogram '" << name << "'" << std::endl;
1702  return nullptr;
1703  }
1704 }
1705 
1706 TH2* AliAnalysisTaskEmcalLight::GetGeneralTH2(const char* name, bool warn)
1707 {
1708  return static_cast<TH2*>(GetGeneralTH1(name, warn));
1709 }
1710 
1711 TProfile* AliAnalysisTaskEmcalLight::GetGeneralTProfile(const char* name, bool warn)
1712 {
1713  return static_cast<TProfile*>(GetGeneralTH1(name, warn));
1714 }
1715 
1717 {
1718  fIsPythia = i;
1719  if (fIsPythia) {
1720  fIsMonteCarlo = kTRUE;
1721  fMCEventHeaderName = "AliGenPythiaEventHeader";
1722  }
1723  else {
1724  if (fMCEventHeaderName == "AliGenPythiaEventHeader") {
1725  fMCEventHeaderName = "";
1726  }
1727  }
1728 }
1729 
1731 {
1732  TClass gen_header_class(name);
1733  if (gen_header_class.InheritsFrom("AliGenEventHeader")) {
1734  fMCEventHeaderName = name;
1735  }
1736  else {
1737  AliWarningStream() << "Class name '" << name << "' does not inherit from 'AliGenEventHeader'. Not setting it." << std::endl;
1738  }
1739 }
Bool_t fSwitchOffLHC15oFaultyBranches
Switch off faulty tree branches in LHC15o AOD trees.
AliClusterContainer * AddClusterContainer(std::string branchName, std::string contName="")
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
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.
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...
TString fCentEst
name of the centrality estimator
UInt_t fTriggerSelectionBitMap
trigger selection bit map
const AliMCParticleIterableContainer all() const
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)
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
AliClusterContainer * GetClusterContainer(std::string name) const
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)
AliParticleContainer * AddParticleContainer(std::string branchName, std::string contName="")
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]
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
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.
AliParticleContainer * GetParticleContainer(std::string name) const
Bool_t fUseXsecFromHeader
!Switch for using cross section from header (if not found in pythia file)
bool Bool_t
Definition: External.C:53
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
Double_t fEPV0C
!event plane V0C
Definition: External.C:196
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)