AliPhysics  06cda10 (06cda10)
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 
18 #include <RVersion.h>
19 #include <TClonesArray.h>
20 #include <TList.h>
21 #include <TObject.h>
22 #include <TH1F.h>
23 #include <TH2F.h>
24 #include <TProfile.h>
25 #include <TSystem.h>
26 #include <TFile.h>
27 #include <TChain.h>
28 #include <TKey.h>
29 
30 #include "AliGenCocktailEventHeader.h"
31 #include "AliStack.h"
32 #include "AliAODEvent.h"
33 #include "AliAnalysisManager.h"
34 #include "AliCentrality.h"
35 #include "AliEMCALGeometry.h"
36 #include "AliESDEvent.h"
37 #include "AliEmcalParticle.h"
38 #include "AliEventplane.h"
39 #include "AliInputEventHandler.h"
40 #include "AliLog.h"
41 #include "AliMCParticle.h"
42 #include "AliVCluster.h"
43 #include "AliVEventHandler.h"
44 #include "AliVParticle.h"
45 #include "AliAODTrack.h"
46 #include "AliVCaloTrigger.h"
47 #include "AliGenPythiaEventHeader.h"
48 #include "AliGenEventHeader.h"
49 #include "AliAODMCHeader.h"
50 #include "AliMCEvent.h"
51 #include "AliEMCALTriggerPatchInfo.h"
52 
53 #include "AliMultSelection.h"
54 
56 
58 
62 
68  fForceBeamType(kNA),
69  fGeneralHistograms(kFALSE),
70  fCreateHisto(kTRUE),
71  fNeedEmcalGeom(kTRUE),
72  fCentBins(),
73  fCentralityEstimation(kNewCentrality),
74  fIsPythia(kFALSE),
75  fIsMonteCarlo(kFALSE),
76  fMCEventHeaderName(),
77  fCaloCellsName(),
78  fCaloTriggersName(),
79  fCaloTriggerPatchInfoName(),
80  fCentEst("V0M"),
81  fParticleCollArray(),
82  fClusterCollArray(),
83  fTriggerSelectionBitMap(0),
84  fMinCent(-1),
85  fMaxCent(-1),
86  fMinVz(-999),
87  fMaxVz(999),
88  fMaxVzDiff(-1),
89  fMinNVertCont(0),
90  fMinPtHard(-1),
91  fMaxPtHard(-1),
92  fMaxMinimumBiasPtHard(-1),
93  fAcceptedTriggerClasses(),
94  fRejectedTriggerClasses(),
95  fMCRejectFilter(kFALSE),
96  fPtHardAndJetPtFactor(0.),
97  fPtHardAndClusterPtFactor(0.),
98  fPtHardAndTrackPtFactor(0.),
99  fSwitchOffLHC15oFaultyBranches(kFALSE),
100  fEventSelectionAfterRun(kFALSE),
101  fSelectGeneratorName(),
102  fMinimumEventWeight(1e-6),
103  fMaximumEventWeight(1e6),
104  fInhibit(kFALSE),
105  fLocalInitialized(kFALSE),
106  fDataType(kAOD),
107  fGeom(0),
108  fCaloCells(0),
109  fCaloTriggers(0),
110  fTriggerPatchInfo(0),
111  fCent(-1),
112  fCentBin(-1),
113  fEPV0(-1.0),
114  fEPV0A(-1.0),
115  fEPV0C(-1.0),
116  fNVertCont(0),
117  fNVertSPDCont(0),
118  fFiredTriggerBitMap(0),
119  fFiredTriggerClasses(),
120  fBeamType(kNA),
121  fMCHeader(0),
122  fPythiaHeader(0),
123  fPtHardBin(0),
124  fPtHard(0),
125  fNTrials(0),
126  fXsection(0),
127  fEventWeight(1),
128  fGeneratorName(),
129  fOutput(0),
130  fHistograms()
131 {
132  fVertex[0] = 0;
133  fVertex[1] = 0;
134  fVertex[2] = 0;
135  fVertexSPD[0] = 0;
136  fVertexSPD[1] = 0;
137  fVertexSPD[2] = 0;
138 }
139 
151  AliAnalysisTaskSE(name),
152  fForceBeamType(kNA),
153  fGeneralHistograms(kFALSE),
154  fCreateHisto(kTRUE),
155  fNeedEmcalGeom(kTRUE),
156  fCentBins(6),
157  fCentralityEstimation(kNewCentrality),
158  fIsPythia(kFALSE),
159  fIsMonteCarlo(kFALSE),
160  fMCEventHeaderName(),
161  fCaloCellsName(),
162  fCaloTriggersName(),
163  fCaloTriggerPatchInfoName(),
164  fCentEst("V0M"),
165  fParticleCollArray(),
166  fClusterCollArray(),
167  fTriggerSelectionBitMap(0),
168  fMinCent(-1),
169  fMaxCent(-1),
170  fMinVz(-999),
171  fMaxVz(999),
172  fMaxVzDiff(-1),
173  fMinNVertCont(0),
174  fMinPtHard(-1),
175  fMaxPtHard(-1),
176  fMaxMinimumBiasPtHard(-1),
177  fAcceptedTriggerClasses(),
178  fRejectedTriggerClasses(),
179  fMCRejectFilter(kFALSE),
180  fPtHardAndJetPtFactor(0.),
181  fPtHardAndClusterPtFactor(0.),
182  fPtHardAndTrackPtFactor(0.),
183  fSwitchOffLHC15oFaultyBranches(kFALSE),
184  fEventSelectionAfterRun(kFALSE),
185  fSelectGeneratorName(),
186  fMinimumEventWeight(1e-6),
187  fMaximumEventWeight(1e6),
188  fInhibit(kFALSE),
189  fLocalInitialized(kFALSE),
190  fDataType(kAOD),
191  fGeom(0),
192  fCaloCells(0),
193  fCaloTriggers(0),
194  fTriggerPatchInfo(0),
195  fCent(0),
196  fCentBin(-1),
197  fEPV0(-1.0),
198  fEPV0A(-1.0),
199  fEPV0C(-1.0),
200  fNVertCont(0),
201  fNVertSPDCont(0),
202  fFiredTriggerBitMap(0),
203  fFiredTriggerClasses(),
204  fBeamType(kNA),
205  fMCHeader(0),
206  fPythiaHeader(0),
207  fPtHardBin(0),
208  fPtHard(0),
209  fNTrials(0),
210  fXsection(0),
211  fEventWeight(1),
212  fGeneratorName(),
213  fOutput(0),
214  fHistograms()
215 {
216  fVertex[0] = 0;
217  fVertex[1] = 0;
218  fVertex[2] = 0;
219  fVertexSPD[0] = 0;
220  fVertexSPD[1] = 0;
221  fVertexSPD[2] = 0;
222 
223  fCentBins[0] = 0;
224  fCentBins[1] = 10;
225  fCentBins[2] = 30;
226  fCentBins[3] = 50;
227  fCentBins[4] = 90;
228  fCentBins[5] = 100;
229 
230  if (fCreateHisto) DefineOutput(1, TList::Class());
231 }
232 
237 {
238  for (auto cont_it : fParticleCollArray) delete cont_it.second;
239  for (auto cont_it : fClusterCollArray) delete cont_it.second;
240 }
241 
262 {
263  if (fInhibit) {
264  AliWarningStream() << "The execution of this task is inhibited. Returning." << std::endl;
265  return;
266  }
267  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
268  if (mgr) {
269  AliVEventHandler *evhand = mgr->GetInputEventHandler();
270  if (evhand) {
271  if (evhand->InheritsFrom("AliESDInputHandler")) {
272  fDataType = kESD;
273  }
274  else {
275  fDataType = kAOD;
276  }
277  }
278  else {
279  AliError("Event handler not found!");
280  }
281  }
282  else {
283  AliError("Analysis manager not found!");
284  }
285 
286  if (!fCreateHisto)
287  return;
288 
289  OpenFile(1);
290  fOutput = new TList();
291  fOutput->SetOwner(); // @suppress("Ambiguous problem")
292 
294 
295  if (!fGeneralHistograms) return;
296 
297  TH1* h = nullptr;
298 
299  if (fIsMonteCarlo) {
300  auto weight_bins = GenerateLogFixedBinArray(1000, fMinimumEventWeight, fMaximumEventWeight, true);
301 
302  h = new TH1F("fHistEventsVsPtHard", "fHistEventsVsPtHard", 1000, 0, 1000);
303  h->GetXaxis()->SetTitle("#it{p}_{T,hard} (GeV/#it{c})");
304  h->GetYaxis()->SetTitle("events");
305  fOutput->Add(h);
306  fHistograms["fHistEventsVsPtHard"] = h;
307 
308  h = new TH1F("fHistTrialsVsPtHard", "fHistTrialsVsPtHard", 1000, 0, 1000);
309  h->GetXaxis()->SetTitle("#it{p}_{T,hard} (GeV/#it{c})");
310  h->GetYaxis()->SetTitle("trials");
311  fOutput->Add(h);
312  fHistograms["fHistTrialsVsPtHard"] = h;
313 
314  h = new TProfile("fHistXsection", "fHistXsection", 50, 0, 50);
315  h->GetXaxis()->SetTitle("#it{p}_{T,hard} bin");
316  h->GetYaxis()->SetTitle("total integrated cross section (mb)");
317  fOutput->Add(h);
318  fHistograms["fHistXsection"] = h;
319 
320  h = new TH1F("fHistXsectionDistribution", "fHistXsectionDistribution", 1000, &weight_bins[0]);
321  h->GetXaxis()->SetTitle("total integrated cross section (mb)");
322  h->GetYaxis()->SetTitle("events");
323  fOutput->Add(h);
324  fHistograms["fHistXsectionDistribution"] = h;
325 
326  h = new TH1F("fHistEventWeights", "fHistEventWeights", 1000, &weight_bins[0]);
327  h->GetXaxis()->SetTitle("weight");
328  h->GetYaxis()->SetTitle("events");
329  fOutput->Add(h);
330  fHistograms["fHistEventWeights"] = h;
331 
332  h = new TH2F("fHistEventWeightsVsPtHard", "fHistEventWeightsVsPtHard", 1000, 0, 1000, 1000, &weight_bins[0]);
333  h->GetXaxis()->SetTitle("#it{p}_{T,hard} (GeV/#it{c})");
334  h->GetYaxis()->SetTitle("event weight");
335  fOutput->Add(h);
336  fHistograms["fHistEventWeightsVsPtHard"] = h;
337 
338  h = new TH1F("fHistEventsVsPtHardNoSel", "fHistEventsVsPtHardNoSel", 1000, 0, 1000);
339  h->GetXaxis()->SetTitle("#it{p}_{T,hard} (GeV/#it{c})");
340  h->GetYaxis()->SetTitle("events");
341  fOutput->Add(h);
342  fHistograms["fHistEventsVsPtHardNoSel"] = h;
343 
344  h = new TH1F("fHistTrialsVsPtHardNoSel", "fHistTrialsVsPtHardNoSel", 1000, 0, 1000);
345  h->GetXaxis()->SetTitle("#it{p}_{T,hard} (GeV/#it{c})");
346  h->GetYaxis()->SetTitle("trials");
347  fOutput->Add(h);
348  fHistograms["fHistTrialsVsPtHardNoSel"] = h;
349 
350  h = new TProfile("fHistXsectionNoSel", "fHistXsectionNoSel", 50, 0, 50);
351  h->GetXaxis()->SetTitle("#it{p}_{T,hard} bin");
352  h->GetYaxis()->SetTitle("total integrated cross section (mb)");
353  fOutput->Add(h);
354  fHistograms["fHistXsectionNoSel"] = h;
355 
356  h = new TH1F("fHistXsectionDistributionNoSel", "fHistXsectionDistributionNoSel", 1000, &weight_bins[0]);
357  h->GetXaxis()->SetTitle("total integrated cross section (mb)");
358  h->GetYaxis()->SetTitle("events");
359  fOutput->Add(h);
360  fHistograms["fHistXsectionDistributionNoSel"] = h;
361 
362  h = new TH1F("fHistEventWeightsNoSel", "fHistEventWeightsNoSel", 1000, &weight_bins[0]);
363  h->GetXaxis()->SetTitle("weight");
364  h->GetYaxis()->SetTitle("events");
365  fOutput->Add(h);
366  fHistograms["fHistEventWeightsNoSel"] = h;
367 
368  h = new TH2F("fHistEventWeightsVsPtHardNoSel", "fHistEventWeightsVsPtHardNoSel", 1000, 0, 1000, 1000, &weight_bins[0]);
369  h->GetXaxis()->SetTitle("#it{p}_{T,hard} (GeV/#it{c})");
370  h->GetYaxis()->SetTitle("event weight");
371  fOutput->Add(h);
372  fHistograms["fHistEventWeightsVsPtHardNoSel"] = h;
373 
374  h = new TH1F("fHistTrialsExternalFile", "fHistTrialsExternalFile", 50, 0, 50);
375  h->GetXaxis()->SetTitle("#it{p}_{T,hard} bin");
376  h->GetYaxis()->SetTitle("trials");
377  fOutput->Add(h);
378  fHistograms["fHistTrialsExternalFile"] = h;
379 
380  h = new TH1F("fHistEventsExternalFile", "fHistEventsExternalFile", 50, 0, 50);
381  h->GetXaxis()->SetTitle("#it{p}_{T,hard} bin");
382  h->GetYaxis()->SetTitle("total events");
383  fOutput->Add(h);
384  fHistograms["fHistEventsExternalFile"] = h;
385 
386  h = new TProfile("fHistXsectionExternalFile", "fHistXsectionExternalFile", 50, 0, 50);
387  h->GetXaxis()->SetTitle("#it{p}_{T,hard} bin");
388  h->GetYaxis()->SetTitle("total integrated cross section (mb)");
389  fOutput->Add(h);
390  fHistograms["fHistXsectionExternalFile"] = h;
391  }
392 
393  h = new TH1F("fHistZVertex","Z vertex position", 60, -30, 30);
394  h->GetXaxis()->SetTitle("V_{#it{z}}");
395  h->GetYaxis()->SetTitle("counts");
396  fOutput->Add(h);
397  fHistograms["fHistZVertex"] = h;
398 
399  h = new TH1F("fHistZVertexNoSel","Z vertex position (no event selection)", 60, -30, 30);
400  h->GetXaxis()->SetTitle("V_{#it{z}}");
401  h->GetYaxis()->SetTitle("counts");
402  fOutput->Add(h);
403  fHistograms["fHistZVertexNoSel"] = h;
404 
406  h = new TH1F("fHistCentrality","Event centrality distribution", 100, 0, 100);
407  h->GetXaxis()->SetTitle("Centrality (%)");
408  h->GetYaxis()->SetTitle("counts");
409  fOutput->Add(h);
410  fHistograms["fHistCentrality"] = h;
411 
412  h = new TH1F("fHistCentralityNoSel","Event centrality distribution (no event selection)", 100, 0, 100);
413  h->GetXaxis()->SetTitle("Centrality (%)");
414  h->GetYaxis()->SetTitle("counts");
415  fOutput->Add(h);
416  fHistograms["fHistCentralityNoSel"] = h;
417  }
418 
419  if (fForceBeamType != kpp) {
420  h = new TH1F("fHistEventPlane","Event plane", 120, -TMath::Pi(), TMath::Pi());
421  h->GetXaxis()->SetTitle("event plane");
422  h->GetYaxis()->SetTitle("counts");
423  fOutput->Add(h);
424  fHistograms["fHistEventPlane"] = h;
425 
426  h = new TH1F("fHistEventPlaneNoSel","Event plane (no event selection)", 120, -TMath::Pi(), TMath::Pi());
427  h->GetXaxis()->SetTitle("event plane");
428  h->GetYaxis()->SetTitle("counts");
429  fOutput->Add(h);
430  fHistograms["fHistEventPlaneNoSel"] = h;
431  }
432 
433  h = new TH1F("fHistEventRejection","Reasons to reject event",30,0,30);
434 #if ROOT_VERSION_CODE < ROOT_VERSION(6,4,2)
435  h->SetBit(TH1::kCanRebin);
436 #else
437  h->SetCanExtend(TH1::kAllAxes);
438 #endif
439  std::array<std::string, 10> labels = {"PhysSel", "Evt Gen Name", "Trg class (acc)", "Trg class (rej)", "Cent", "vertex contr.", "Vz", "VzSPD", "SelPtHardBin", "MCOutlier"};
440  int i = 1;
441  for (auto label : labels) {
442  h->GetXaxis()->SetBinLabel(i, label.c_str());
443  i++;
444  }
445  h->GetYaxis()->SetTitle("counts");
446  fOutput->Add(h);
447  fHistograms["fHistEventRejection"] = h;
448 
449  h = new TH1F("fHistTriggerClasses","fHistTriggerClasses",3,0,3);
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  fOutput->Add(h);
456  fHistograms["fHistTriggerClasses"] = h;
457 
458  h = new TH1F("fHistTriggerClassesNoSel","fHistTriggerClassesNoSel",3,0,3);
459 #if ROOT_VERSION_CODE < ROOT_VERSION(6,4,2)
460  h->SetBit(TH1::kCanRebin);
461 #else
462  h->SetCanExtend(TH1::kAllAxes);
463 #endif
464  fOutput->Add(h);
465  fHistograms["fHistTriggerClassesNoSel"] = h;
466 
467  h = new TH1F("fHistEventCount","fHistEventCount",2,0,2);
468  h->GetXaxis()->SetBinLabel(1,"Accepted");
469  h->GetXaxis()->SetBinLabel(2,"Rejected");
470  h->GetYaxis()->SetTitle("counts");
471  fOutput->Add(h);
472  fHistograms["fHistEventCount"] = h;
473 
474  PostData(1, fOutput);
475 }
476 
492 {
493  if (eventSelected) {
494  if (fIsMonteCarlo) {
495  GetGeneralTH1("fHistEventsVsPtHard", true)->Fill(fPtHard);
496  GetGeneralTH1("fHistTrialsVsPtHard", true)->Fill(fPtHard, fNTrials);
497  GetGeneralTH1("fHistEventWeights", true)->Fill(fEventWeight);
498  GetGeneralTH2("fHistEventWeightsVsPtHard", true)->Fill(fPtHard, fEventWeight);
499  GetGeneralTH1("fHistXsectionDistribution", true)->Fill(fXsection);
500  GetGeneralTProfile("fHistXsection", true)->Fill(fPtHardBin, fXsection);
501  }
502 
503  GetGeneralTH1("fHistZVertex")->Fill(fVertex[2]);
504 
505  TH1* hCent = GetGeneralTH1("fHistCentrality");
506  if (hCent) hCent->Fill(fCent);
507 
508  TH1* hEventPlane = GetGeneralTH1("fHistEventPlane");
509  if (hEventPlane) hEventPlane->Fill(fEPV0);
510 
511  TH1* hTriggerClasses = GetGeneralTH1("fHistTriggerClasses");
512  for (auto fired_trg : fFiredTriggerClasses) hTriggerClasses->Fill(fired_trg.c_str(), 1);
513  }
514  else {
515  if (fIsMonteCarlo) {
516  GetGeneralTH1("fHistEventsVsPtHardNoSel", true)->Fill(fPtHard);
517  GetGeneralTH1("fHistTrialsVsPtHardNoSel", true)->Fill(fPtHard, fNTrials);
518  GetGeneralTH1("fHistEventWeightsNoSel", true)->Fill(fEventWeight);
519  GetGeneralTH2("fHistEventWeightsVsPtHardNoSel", true)->Fill(fPtHard, fEventWeight);
520  GetGeneralTH1("fHistXsectionDistributionNoSel", true)->Fill(fXsection);
521  GetGeneralTProfile("fHistXsectionNoSel", true)->Fill(fPtHardBin, fXsection);
522  }
523 
524  GetGeneralTH1("fHistZVertexNoSel", true)->Fill(fVertex[2]);
525 
526  TH1* hCent = GetGeneralTH1("fHistCentralityNoSel");
527  if (hCent) hCent->Fill(fCent);
528 
529  TH1* hEventPlane = GetGeneralTH1("fHistEventPlaneNoSel");
530  if (hEventPlane) hEventPlane->Fill(fEPV0);
531 
532  TH1* hTriggerClasses = GetGeneralTH1("fHistTriggerClassesNoSel", true);
533  for (auto fired_trg : fFiredTriggerClasses) hTriggerClasses->Fill(fired_trg.c_str(), 1);
534  }
535 
536  return kTRUE;
537 }
538 
559 {
560  if (fInhibit) {
561  AliWarningStream() << "The execution of this task is inhibited. Returning." << std::endl;
562  return;
563  }
564 
565  if (!fLocalInitialized) ExecOnce();
566 
567  if (!fLocalInitialized) return;
568 
569  if (!RetrieveEventObjects()) return;
570 
571  Bool_t eventSelected = IsEventSelected();
572 
574  if (eventSelected) {
575  GetGeneralTH1("fHistEventCount", true)->Fill("Accepted",1);
576  }
577  else {
578  GetGeneralTH1("fHistEventCount", true)->Fill("Rejected",1);
579  }
580 
581  FillGeneralHistograms(kFALSE);
582  if (eventSelected) FillGeneralHistograms(kTRUE);
583  }
584 
585  Bool_t runOk = kFALSE;
586  if (eventSelected || fEventSelectionAfterRun) runOk = Run();
587 
588  if (fCreateHisto && eventSelected && runOk) FillHistograms();
589 
590  if (fCreateHisto && fOutput) {
591  // information for this iteration of the UserExec in the container
592  PostData(1, fOutput);
593  }
594 }
595 
607 Bool_t AliAnalysisTaskEmcalLight::PythiaInfoFromFile(const char* currFile, Float_t &xsec, Float_t &trials, Int_t &pthard)
608 {
609 
610  TString file(currFile);
611  xsec = 0;
612  trials = 1;
613 
614  if (file.Contains(".zip#")) {
615  Ssiz_t pos1 = file.Index("root_archive",12,0,TString::kExact);
616  Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
617  Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
618  file.Replace(pos+1,pos2-pos1,"");
619  } else {
620  // not an archive take the basename....
621  file.ReplaceAll(gSystem->BaseName(file.Data()),"");
622  }
623  AliDebug(1,Form("File name: %s",file.Data()));
624 
625  // Get the pt hard bin
626  TString strPthard(file);
627 
628  strPthard.Remove(strPthard.Last('/'));
629  strPthard.Remove(strPthard.Last('/'));
630  if (strPthard.Contains("AOD")) strPthard.Remove(strPthard.Last('/'));
631  strPthard.Remove(0,strPthard.Last('/')+1);
632  if (strPthard.IsDec()) {
633  pthard = strPthard.Atoi();
634  }
635  else {
636  AliWarning(Form("Could not extract file number from path %s", strPthard.Data()));
637  pthard = -1;
638  }
639 
640  // 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
641  TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root"));
642 
643  if (!fxsec) {
644  // next trial fetch the histgram file
645  fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
646  if (!fxsec) {
647  // not a severe condition but inciate that we have no information
648  return kFALSE;
649  } else {
650  // find the tlist we want to be independtent of the name so use the Tkey
651  TKey* key = static_cast<TKey*>(fxsec->GetListOfKeys()->At(0));
652  if (!key) {
653  fxsec->Close();
654  return kFALSE;
655  }
656  TList *list = dynamic_cast<TList*>(key->ReadObj());
657  if (!list) {
658  fxsec->Close();
659  return kFALSE;
660  }
661  xsec = static_cast<TProfile*>(list->FindObject("h1Xsec"))->GetBinContent(1);
662  trials = static_cast<TH1F*>(list->FindObject("h1Trials"))->GetBinContent(1);
663  fxsec->Close();
664  }
665  } else { // no tree pyxsec.root
666  TTree *xtree = static_cast<TTree*>(fxsec->Get("Xsection"));
667  if (!xtree) {
668  fxsec->Close();
669  return kFALSE;
670  }
671  UInt_t ntrials = 0;
672  Double_t xsection = 0;
673  xtree->SetBranchAddress("xsection",&xsection);
674  xtree->SetBranchAddress("ntrials",&ntrials);
675  xtree->GetEntry(0);
676  trials = ntrials;
677  xsec = xsection;
678  fxsec->Close();
679  }
680  return kTRUE;
681 }
682 
697 {
698  if (!fIsPythia) return kTRUE;
699 
700  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
701  if (!tree) {
702  AliError(Form("%s - UserNotify: No current tree!",GetName()));
703  return kFALSE;
704  }
705 
706  Float_t xsection = 0;
707  Float_t trials = 0;
708  Int_t pthardbin = 0;
709 
710  TFile *curfile = tree->GetCurrentFile();
711  if (!curfile) {
712  AliError(Form("%s - UserNotify: No current file!",GetName()));
713  return kFALSE;
714  }
715 
716  TChain *chain = dynamic_cast<TChain*>(tree);
717  if (chain) tree = chain->GetTree();
718 
719  Int_t nevents = tree->GetEntriesFast();
720 
721  Bool_t res = PythiaInfoFromFile(curfile->GetName(), xsection, trials, pthardbin);
722 
723  fPtHardBin = pthardbin >= 0 ? pthardbin : 0;
724 
725  if (!res) return kTRUE;
726 
728  GetGeneralTH1("fHistTrialsExternalFile", true)->Fill(fPtHardBin, trials);
729  GetGeneralTProfile("fHistXsectionExternalFile", true)->Fill(fPtHardBin, xsection);
730  GetGeneralTH1("fHistEventsExternalFile", true)->Fill(fPtHardBin, nevents);
731  }
732 
733  return kTRUE;
734 }
735 
747 {
748  if (!InputEvent()) {
749  AliError(Form("%s: Could not retrieve event! Returning!", GetName()));
750  return;
751  }
752 
753  if (fNeedEmcalGeom && !fGeom) {
754  fGeom = AliEMCALGeometry::GetInstanceFromRunNumber(InputEvent()->GetRunNumber());
755  if (!fGeom) {
756  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()));
757  return;
758  }
759  }
760 
761  if (fSwitchOffLHC15oFaultyBranches && dynamic_cast<AliAODEvent*>(InputEvent())) {
762  TTree *aodTree = AliAnalysisManager::GetAnalysisManager()->GetTree();
763  aodTree->SetBranchStatus("D0toKpi.fPx", 0);
764  aodTree->SetBranchStatus("D0toKpi.fPy", 0);
765  aodTree->SetBranchStatus("D0toKpi.fPz", 0);
766  aodTree->SetBranchStatus("D0toKpi.fd0", 0);
767  aodTree->SetBranchStatus("Charm3Prong.fPx", 0);
768  aodTree->SetBranchStatus("Charm3Prong.fPy", 0);
769  aodTree->SetBranchStatus("Charm3Prong.fPz", 0);
770  aodTree->SetBranchStatus("Charm3Prong.fd0", 0);
771  aodTree->SetBranchStatus("Dstar.fPx", 0);
772  aodTree->SetBranchStatus("Dstar.fPy", 0);
773  aodTree->SetBranchStatus("Dstar.fPz", 0);
774  aodTree->SetBranchStatus("Dstar.fd0", 0);
775  }
776 
777  //Load all requested track branches - each container knows name already
778  for (auto cont_it : fParticleCollArray) cont_it.second->SetArray(InputEvent());
779 
780  //Load all requested cluster branches - each container knows name already
781  for (auto cont_it : fClusterCollArray) cont_it.second->SetArray(InputEvent());
782 
783  if (!fCaloCellsName.IsNull() && !fCaloCells) {
784  fCaloCells = dynamic_cast<AliVCaloCells*>(InputEvent()->FindListObject(fCaloCellsName));
785  if (!fCaloCells) {
786  AliError(Form("%s: Could not retrieve cells %s!", GetName(), fCaloCellsName.Data()));
787  return;
788  }
789  }
790 
791  if (!fCaloTriggersName.IsNull() && !fCaloTriggers) {
792  fCaloTriggers = dynamic_cast<AliVCaloTrigger*>(InputEvent()->FindListObject(fCaloTriggersName));
793  if (!fCaloTriggers) {
794  AliError(Form("%s: Could not retrieve calo triggers %s!", GetName(), fCaloTriggersName.Data()));
795  return;
796  }
797  }
798 
799  if (!fCaloTriggerPatchInfoName.IsNull() && !fTriggerPatchInfo) {
800  fTriggerPatchInfo = GetArrayFromEvent(fCaloTriggerPatchInfoName.Data(),"AliEMCALTriggerPatchInfo");
801  if (!fTriggerPatchInfo) {
802  AliError(Form("%s: Could not retrieve calo trigger patch info %s!", GetName(), fCaloTriggerPatchInfoName.Data()));
803  return;
804  }
805 
806  }
807 
808  fLocalInitialized = kTRUE;
809 }
810 
817 {
818  if (fForceBeamType != kNA)
819  return fForceBeamType;
820 
821  AliESDEvent *esd = dynamic_cast<AliESDEvent*>(InputEvent());
822  if (esd) {
823  const AliESDRun *run = esd->GetESDRun();
824  TString beamType = run->GetBeamType();
825  if (beamType == "p-p")
826  return kpp;
827  else if (beamType == "A-A")
828  return kAA;
829  else if (beamType == "p-A")
830  return kpA;
831  else
832  return kNA;
833  } else {
834  Int_t runNumber = InputEvent()->GetRunNumber();
835  // All run number ranges taken from the RCT
836  if ((runNumber >= 136833 && runNumber <= 139517) || // LHC10h
837  (runNumber >= 167693 && runNumber <= 170593) || // LHC11h
838  (runNumber >= 244824 && runNumber <= 246994)) { // LHC15o
839  return kAA;
840  } else if ((runNumber >= 188356 && runNumber <= 188366) || // LHC12g
841  (runNumber >= 195164 && runNumber <= 197388) || // LHC13b-f
842  (runNumber >= 265015 && runNumber <= 267166)) { // LHC16q-t
843  return kpA;
844  } else {
845  return kpp;
846  }
847  }
848 }
849 
872 {
873  TH1* hEventRejection = GetGeneralTH1("fHistEventRejection", true);
874 
876  if (fGeneralHistograms) hEventRejection->Fill("PhysSel",1);
877  return kFALSE;
878  }
879 
880  if (!fSelectGeneratorName.IsNull() && !fGeneratorName.IsNull()) {
881  if (!fGeneratorName.Contains(fSelectGeneratorName)) {
882  if (fGeneralHistograms) hEventRejection->Fill("Evt Gen Name",1);
883  return kFALSE;
884  }
885  }
886 
887  Bool_t acceptedTrgClassFound = kFALSE;
888  if (fAcceptedTriggerClasses.size() > 0) {
889  for (auto acc_trg : fAcceptedTriggerClasses) {
890  std::string teststring(acc_trg);
891  bool fullmatch(false);
892  auto posexact = acc_trg.find("EXACT");
893  if(posexact != std::string::npos) {
894  fullmatch = true;
895  teststring.erase(posexact, 5);
896  }
897  for (auto fired_trg : fFiredTriggerClasses) {
898  bool classmatch = fullmatch ? teststring == fired_trg : fired_trg.find(teststring) != std::string::npos;
899  if (classmatch) {
900  acceptedTrgClassFound = kTRUE;
901  break;
902  }
903  }
904  if (acceptedTrgClassFound) break;
905  }
906 
907  if (!acceptedTrgClassFound) {
908  if (fGeneralHistograms) hEventRejection->Fill("Trg class (acc)",1);
909  return kFALSE;
910  }
911  }
912 
913  if (fRejectedTriggerClasses.size() > 0) {
914  for (auto rej_trg : fRejectedTriggerClasses) {
915  std::string teststring(rej_trg);
916  bool fullmatch(false);
917  auto posexact = rej_trg.find("EXACT");
918  if(posexact != std::string::npos) {
919  fullmatch = true;
920  teststring.erase(posexact, 5);
921  }
922  for (auto fired_trg : fFiredTriggerClasses) {
923  bool classmatch = fullmatch ? teststring == fired_trg : fired_trg.find(teststring) != std::string::npos;
924  if (classmatch) {
925  if (fGeneralHistograms) hEventRejection->Fill("Trg class (rej)",1);
926  return kFALSE;
927  }
928  }
929  }
930  }
931 
932  if (fMinCent < fMaxCent && fMaxCent > 0) {
933  if (fCent < fMinCent || fCent > fMaxCent) {
934  if (fGeneralHistograms) hEventRejection->Fill("Cent",1);
935  return kFALSE;
936  }
937  }
938 
939  if (fNVertCont < fMinNVertCont) {
940  if (fGeneralHistograms) hEventRejection->Fill("vertex contr.",1);
941  return kFALSE;
942  }
943 
944  if (fMinVz < fMaxVz) {
945  if (fVertex[2] < fMinVz || fVertex[2] > fMaxVz) {
946  if (fGeneralHistograms) hEventRejection->Fill("Vz",1);
947  return kFALSE;
948  }
949  }
950 
951  if (fMaxVzDiff >= 0) {
952  if (fNVertSPDCont > 0) {
953  Double_t vzSPD = fVertexSPD[2];
954  Double_t dvertex = TMath::Abs(fVertex[2] - vzSPD);
955  //if difference larger than fZvertexDiff
956  if (dvertex > fMaxVzDiff) {
957  if (fGeneralHistograms) hEventRejection->Fill("VzSPD",1);
958  return kFALSE;
959  }
960  }
961  }
962 
963  if (fMinPtHard >= 0 && fPtHard < fMinPtHard) {
964  if (fGeneralHistograms) hEventRejection->Fill("SelPtHardBin",1);
965  return kFALSE;
966  }
967 
968  if (fMaxPtHard >= 0 && fPtHard >= fMaxPtHard) {
969  if (fGeneralHistograms) hEventRejection->Fill("SelPtHardBin",1);
970  return kFALSE;
971  }
972 
974  if (fGeneralHistograms) hEventRejection->Fill("SelPtHardBin",1);
975  return kFALSE;
976  }
977 
978  // Reject filter for MC data
979  if (!CheckMCOutliers()) {
980  if (fGeneralHistograms) hEventRejection->Fill("MCOutlier",1);
981  return kFALSE;
982  }
983 
984  return kTRUE;
985 }
986 
995 TClonesArray *AliAnalysisTaskEmcalLight::GetArrayFromEvent(const char *name, const char *clname)
996 {
997  TClonesArray *arr = 0;
998  TString sname(name);
999  if (!sname.IsNull()) {
1000  arr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(sname));
1001  if (!arr) {
1002  AliWarning(Form("%s: Could not retrieve array with name %s!", GetName(), name));
1003  return 0;
1004  }
1005  } else {
1006  return 0;
1007  }
1008 
1009  if (!clname)
1010  return arr;
1011 
1012  TString objname(arr->GetClass()->GetName());
1013  TClass cls(objname);
1014  if (!cls.InheritsFrom(clname)) {
1015  AliWarning(Form("%s: Objects of type %s in %s are not inherited from %s!",
1016  GetName(), cls.GetName(), name, clname));
1017  return 0;
1018  }
1019  return arr;
1020 }
1021 
1027 {
1028  fVertex[0] = 0;
1029  fVertex[1] = 0;
1030  fVertex[2] = 0;
1031  fNVertCont = 0;
1032 
1033  fVertexSPD[0] = 0;
1034  fVertexSPD[1] = 0;
1035  fVertexSPD[2] = 0;
1036  fNVertSPDCont = 0;
1037 
1038  fFiredTriggerClasses.clear();
1039  std::stringstream firedClasses(InputEvent()->GetFiredTriggerClasses().Data());
1040  while (firedClasses.good()) {
1041  std::string trgClass;
1042  firedClasses >> trgClass;
1043  if (!trgClass.empty()) fFiredTriggerClasses.push_back(trgClass);
1044  }
1045 
1046  if (fDataType == kESD) {
1047  fFiredTriggerBitMap = static_cast<AliInputEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler())->IsEventSelected();
1048  }
1049  else {
1050  fFiredTriggerBitMap = static_cast<AliVAODHeader*>(InputEvent()->GetHeader())->GetOfflineTrigger();
1051  }
1052 
1053  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
1054  if (vert) {
1055  vert->GetXYZ(fVertex);
1056  fNVertCont = vert->GetNContributors();
1057  }
1058 
1059  const AliVVertex *vertSPD = InputEvent()->GetPrimaryVertexSPD();
1060  if (vertSPD) {
1061  vertSPD->GetXYZ(fVertexSPD);
1062  fNVertSPDCont = vertSPD->GetNContributors();
1063  }
1064 
1065  fBeamType = GetBeamType();
1066 
1067  fCent = 99;
1068  fCentBin = -1;
1069  fEPV0 = -999;
1070  fEPV0A = -999;
1071  fEPV0C = -999;
1072 
1074  // New centrality estimation (AliMultSelection)
1075  // See https://twiki.cern.ch/twiki/bin/viewauth/ALICE/AliMultSelectionCalibStatus for calibration status period-by-period)
1076  AliMultSelection *MultSelection = static_cast<AliMultSelection*>(InputEvent()->FindListObject("MultSelection"));
1077  if (MultSelection) {
1078  fCent = MultSelection->GetMultiplicityPercentile(fCentEst.Data());
1079  }
1080  else {
1081  AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
1082  }
1083  }
1084  else if (fCentralityEstimation == kOldCentrality) {
1085  // Old centrality estimation (AliCentrality, works only on Run-1 PbPb and pPb)
1086  AliCentrality *aliCent = InputEvent()->GetCentrality();
1087  if (aliCent) {
1088  fCent = aliCent->GetCentralityPercentile(fCentEst.Data());
1089  }
1090  else {
1091  AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
1092  }
1093  }
1094  if (!fCentBins.empty() && fCentralityEstimation != kNoCentrality) {
1095  for (auto cent_it = fCentBins.begin(); cent_it != fCentBins.end() - 1; cent_it++) {
1096  if (fCent >= *cent_it && fCent < *(cent_it+1)) fCentBin = cent_it - fCentBins.begin();
1097  }
1098  }
1099  else {
1100  fCentBin = 0;
1101  }
1102 
1103  if (fBeamType == kAA || fBeamType == kpA ) {
1104  AliEventplane *aliEP = InputEvent()->GetEventplane();
1105  if (aliEP) {
1106  fEPV0 = aliEP->GetEventplane("V0" ,InputEvent());
1107  fEPV0A = aliEP->GetEventplane("V0A",InputEvent());
1108  fEPV0C = aliEP->GetEventplane("V0C",InputEvent());
1109  } else {
1110  AliWarning(Form("%s: Could not retrieve event plane information!", GetName()));
1111  }
1112  }
1113 
1114  if (fIsMonteCarlo && MCEvent()) {
1115  AliGenEventHeader* header = MCEvent()->GenEventHeader();
1116  if (fMCEventHeaderName.IsNull()) {
1117  fMCHeader = header;
1118  }
1119  else {
1120  if (header->InheritsFrom(fMCEventHeaderName)) {
1121  fMCHeader = header;
1122  }
1123  else if (header->InheritsFrom("AliGenCocktailEventHeader")) {
1124  AliGenCocktailEventHeader* cocktailHeader = static_cast<AliGenCocktailEventHeader*>(header);
1125  TList* headers = cocktailHeader->GetHeaders();
1126  for (auto obj : *headers) { // @suppress("Symbol is not resolved")
1127  if (obj->InheritsFrom(fMCEventHeaderName)){
1128  fMCHeader = static_cast<AliGenEventHeader*>(obj);
1129  break;
1130  }
1131  }
1132  }
1133  }
1134  if (fMCHeader) {
1135  fEventWeight = fMCHeader->EventWeight();
1136  if (fIsPythia) {
1137  fPythiaHeader = static_cast<AliGenPythiaEventHeader*>(fMCHeader);
1138  fPtHard = fPythiaHeader->GetPtHard();
1139  fXsection = fPythiaHeader->GetXsection();
1140  fNTrials = fPythiaHeader->Trials();
1141  }
1142  }
1143  }
1144 
1145  for (auto cont_it : fParticleCollArray) cont_it.second->NextEvent(InputEvent());
1146  for (auto cont_it : fClusterCollArray) cont_it.second->NextEvent(InputEvent());
1147 
1148  return kTRUE;
1149 }
1150 
1159 AliParticleContainer* AliAnalysisTaskEmcalLight::AddParticleContainer(std::string branchName, std::string contName)
1160 {
1161  if (branchName.size() == 0) return 0;
1162 
1163  AliParticleContainer* cont = 0;
1164 
1165  if (branchName == "tracks" || branchName == "Tracks") cont = new AliTrackContainer(branchName.c_str());
1166  else if (branchName == "mcparticles") cont = new AliMCParticleContainer(branchName.c_str());
1167  else cont = new AliParticleContainer(branchName.c_str());
1168 
1169  if (contName.size() > 0) cont->SetName(contName.c_str());
1170 
1171  AdoptParticleContainer(cont);
1172 
1173  return cont;
1174 }
1175 
1184 AliClusterContainer* AliAnalysisTaskEmcalLight::AddClusterContainer(std::string branchName, std::string contName)
1185 {
1186  if (branchName.size() == 0) return 0;
1187 
1188  AliClusterContainer* cont = new AliClusterContainer(branchName.c_str());
1189 
1190  if (contName.size() > 0) cont->SetName(contName.c_str());
1191 
1192  AdoptClusterContainer(cont);
1193 
1194  return cont;
1195 }
1196 
1203 {
1204  std::map<std::string, AliParticleContainer*>::const_iterator cont_it = fParticleCollArray.find(name);
1205  if (cont_it != fParticleCollArray.end()) return cont_it->second;
1206  else return nullptr;
1207 }
1208 
1215 {
1216  std::map<std::string, AliClusterContainer*>::const_iterator cont_it = fClusterCollArray.find(name);
1217  if (cont_it != fClusterCollArray.end()) return cont_it->second;
1218  else return nullptr;
1219 }
1220 
1227 {
1228  if (!(InputEvent()->FindListObject(obj->GetName()))) {
1229  InputEvent()->AddObject(obj);
1230  }
1231  else {
1232  if (!attempt) {
1233  AliFatal(Form("%s: Container with name %s already present. Aborting", GetName(), obj->GetName()));
1234  }
1235  }
1236 }
1237 
1246 {
1247 
1248  if (!fGeom) {
1249  AliWarning(Form("%s - AliAnalysisTaskEmcalBase::IsTrackInEmcalAcceptance - Geometry is not available!", GetName()));
1250  return kFALSE;
1251  }
1252 
1253  Double_t minPhi = fGeom->GetArm1PhiMin() - edges;
1254  Double_t maxPhi = fGeom->GetArm1PhiMax() + edges;
1255 
1256  if (part->Phi() > minPhi && part->Phi() < maxPhi) {
1257  return kTRUE;
1258  }
1259  else {
1260  return kFALSE;
1261  }
1262 }
1263 
1265 {
1266  axis->SetBinLabel(1, "NullObject");
1267  axis->SetBinLabel(2, "Pt");
1268  axis->SetBinLabel(3, "Acceptance");
1269  axis->SetBinLabel(4, "MCLabel");
1270  axis->SetBinLabel(5, "BitMap");
1271  axis->SetBinLabel(6, "HF cut");
1272  axis->SetBinLabel(7, "Bit6");
1273  axis->SetBinLabel(8, "NotHybridTrack");
1274  axis->SetBinLabel(9, "MCFlag");
1275  axis->SetBinLabel(10, "MCGenerator");
1276  axis->SetBinLabel(11, "ChargeCut");
1277  axis->SetBinLabel(12, "MinDistanceTPCSectorEdge");
1278  axis->SetBinLabel(13, "Bit12");
1279  axis->SetBinLabel(14, "IsEMCal");
1280  axis->SetBinLabel(15, "Time");
1281  axis->SetBinLabel(16, "Energy");
1282  axis->SetBinLabel(17, "ExoticCut");
1283  axis->SetBinLabel(18, "Bit17");
1284  axis->SetBinLabel(19, "Area");
1285  axis->SetBinLabel(20, "AreaEmc");
1286  axis->SetBinLabel(21, "ZLeadingCh");
1287  axis->SetBinLabel(22, "ZLeadingEmc");
1288  axis->SetBinLabel(23, "NEF");
1289  axis->SetBinLabel(24, "MinLeadPt");
1290  axis->SetBinLabel(25, "MaxTrackPt");
1291  axis->SetBinLabel(26, "MaxClusterPt");
1292  axis->SetBinLabel(27, "Flavour");
1293  axis->SetBinLabel(28, "TagStatus");
1294  axis->SetBinLabel(29, "MinNConstituents");
1295  axis->SetBinLabel(30, "Bit29");
1296  axis->SetBinLabel(31, "Bit30");
1297  axis->SetBinLabel(32, "Bit31");
1298 }
1299 
1306 Double_t AliAnalysisTaskEmcalLight::GetParallelFraction(AliVParticle* part1, AliVParticle* part2)
1307 {
1308  TVector3 vect1(part1->Px(), part1->Py(), part1->Pz());
1309  TVector3 vect2(part2->Px(), part2->Py(), part2->Pz());
1310  Double_t z = (vect1 * vect2) / (vect2 * vect2);
1311  return z;
1312 }
1313 
1320 Double_t AliAnalysisTaskEmcalLight::GetParallelFraction(const TVector3& vect1, AliVParticle* part2)
1321 {
1322  TVector3 vect2(part2->Px(), part2->Py(), part2->Pz());
1323  Double_t z = (vect1 * vect2) / (vect2 * vect2);
1324  return z;
1325 }
1326 
1335 void AliAnalysisTaskEmcalLight::GetEtaPhiDiff(const AliVTrack *t, const AliVCluster *v, Double_t &phidiff, Double_t &etadiff)
1336 {
1337  phidiff = 999;
1338  etadiff = 999;
1339 
1340  if (!t||!v) return;
1341 
1342  Double_t veta = t->GetTrackEtaOnEMCal();
1343  Double_t vphi = t->GetTrackPhiOnEMCal();
1344 
1345  Float_t pos[3] = {0};
1346  v->GetPosition(pos);
1347  TVector3 cpos(pos);
1348  Double_t ceta = cpos.Eta();
1349  Double_t cphi = cpos.Phi();
1350  etadiff=veta-ceta;
1351  phidiff=TVector2::Phi_mpi_pi(vphi-cphi);
1352 }
1353 
1360 {
1361  Byte_t ret = 0;
1362  if (t->TestBit(BIT(22)) && !t->TestBit(BIT(23)))
1363  ret = 1;
1364  else if (!t->TestBit(BIT(22)) && t->TestBit(BIT(23)))
1365  ret = 2;
1366  else if (t->TestBit(BIT(22)) && t->TestBit(BIT(23)))
1367  ret = 3;
1368  return ret;
1369 }
1370 
1380 Byte_t AliAnalysisTaskEmcalLight::GetTrackType(const AliAODTrack *aodTrack, UInt_t filterBit1, UInt_t filterBit2)
1381 {
1382 
1383  Int_t res = 0;
1384 
1385  if (aodTrack->TestFilterBit(filterBit1)) {
1386  res = 0;
1387  }
1388  else if (aodTrack->TestFilterBit(filterBit2)) {
1389  if ((aodTrack->GetStatus()&AliVTrack::kITSrefit)!=0) {
1390  res = 1;
1391  }
1392  else {
1393  res = 2;
1394  }
1395  }
1396  else {
1397  res = 3;
1398  }
1399 
1400  return res;
1401 }
1402 
1409 {
1410  EBeamType_t b = kpp;
1411  if ((runnumber >= 136833 && runnumber <= 139517) || // LHC10h Run-1 (Pb-Pb)
1412  (runnumber >= 167693 && runnumber <= 170593) || // LHC11h Run-1 (Pb-Pb)
1413  (runnumber >= 244824 && runnumber <= 246994)) { // LHC15o Run-2 (Pb-Pb)
1414  b = kAA;
1415  }
1416  else if ((runnumber > 188356 && runnumber <= 188503) || // LHC12g Run-1 (p-Pb pilot)
1417  (runnumber >= 195164 && runnumber <= 197388) || // LHC13b,c,d,e,f Run-1 (p-Pb)
1418  (runnumber >= 265077 && runnumber <= 267166)) { // LHC16 Run-2 (p-Pb)
1419  b = kpA;
1420  }
1421  return b;
1422 }
1423 
1430 {
1431  if (!fPythiaHeader || !fMCRejectFilter) return kTRUE;
1432 
1433  // Condition 1: Pythia jet / pT-hard > factor
1434  if (fPtHardAndJetPtFactor > 0.) {
1435  AliTLorentzVector jet;
1436 
1437  Int_t nTriggerJets = fPythiaHeader->NTriggerJets();
1438 
1439  AliDebug(1,Form("Njets: %d, pT Hard %f",nTriggerJets, fPtHard));
1440 
1441  Float_t tmpjet[]={0,0,0,0};
1442  for (Int_t ijet = 0; ijet< nTriggerJets; ijet++) {
1443  fPythiaHeader->TriggerJet(ijet, tmpjet);
1444 
1445  jet.SetPxPyPzE(tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3]);
1446 
1447  AliDebug(1,Form("jet %d; pycell jet pT %f",ijet, jet.Pt()));
1448 
1449  //Compare jet pT and pt Hard
1450  if (jet.Pt() > fPtHardAndJetPtFactor * fPtHard) {
1451  AliInfo(Form("Reject jet event with : pT Hard %2.2f, pycell jet pT %2.2f, rejection factor %1.1f\n", fPtHard, jet.Pt(), fPtHardAndJetPtFactor));
1452  return kFALSE;
1453  }
1454  }
1455  }
1456  // end condition 1
1457 
1458  // Condition 2 : Reconstructed EMCal cluster pT / pT-hard > factor
1459  if (fPtHardAndClusterPtFactor > 0.) {
1460  AliClusterContainer* mccluscont = GetClusterContainer(0);
1461  if ((Bool_t)mccluscont) {
1462  for (auto cluster : mccluscont->all()) {// Not cuts applied ; use accept for cuts
1463  Float_t ecluster = cluster->E();
1464 
1465  if (ecluster > (fPtHardAndClusterPtFactor * fPtHard)) {
1466  AliInfo(Form("Reject : ecluster %2.2f, calo %d, factor %2.2f, ptHard %f",ecluster,cluster->GetType(),fPtHardAndClusterPtFactor,fPtHard));
1467  return kFALSE;
1468  }
1469  }
1470  }
1471  }
1472  // end condition 2
1473 
1474  // condition 3 : Reconstructed track pT / pT-hard >factor
1475  if (fPtHardAndTrackPtFactor > 0.) {
1476  AliMCParticleContainer* mcpartcont = dynamic_cast<AliMCParticleContainer*>(GetParticleContainer(0));
1477  if ((Bool_t)mcpartcont) {
1478  for (auto mctrack : mcpartcont->all()) {// Not cuts applied ; use accept for cuts
1479  Float_t trackpt = mctrack->Pt();
1480  if (trackpt > (fPtHardAndTrackPtFactor * fPtHard) ) {
1481  AliInfo(Form("Reject : track %2.2f, factor %2.2f, ptHard %f", trackpt, fPtHardAndTrackPtFactor, fPtHard));
1482  return kFALSE;
1483  }
1484  }
1485  }
1486  }
1487  // end condition 3
1488 
1489  return kTRUE;
1490 }
1491 
1501 {
1502  Double_t dphi = -999;
1503  const Double_t tpi = TMath::TwoPi();
1504 
1505  if (phia < 0) phia += tpi;
1506  else if (phia > tpi) phia -= tpi;
1507  if (phib < 0) phib += tpi;
1508  else if (phib > tpi) phib -= tpi;
1509  dphi = phib - phia;
1510  if (dphi < rangeMin) dphi += tpi;
1511  else if (dphi > rangeMax) dphi -= tpi;
1512 
1513  return dphi;
1514 }
1515 
1525 void AliAnalysisTaskEmcalLight::GenerateFixedBinArray(int n, double min, double max, std::vector<double>& array, bool last)
1526 {
1527  double binWidth = (max - min) / n;
1528  double v = min;
1529  if (last) n++;
1530  for (int i = 0; i < n; i++) {
1531  array.push_back(v);
1532  v += binWidth;
1533  }
1534 }
1535 
1545 std::vector<double> AliAnalysisTaskEmcalLight::GenerateFixedBinArray(int n, double min, double max, bool last)
1546 {
1547  std::vector<double> array;
1548  GenerateFixedBinArray(n, min, max, array, last);
1549  return array;
1550 }
1551 
1561 void AliAnalysisTaskEmcalLight::GenerateLogFixedBinArray(int n, double min, double max, std::vector<double>& array, bool last)
1562 {
1563  if (min <= 0 || max < min) {
1564  AliErrorClassStream() << "Cannot generate a log scale fixed-bin array with limits " << min << ", " << max << std::endl;
1565  return;
1566  }
1567  double binWidth = std::pow(max / min, 1.0 / n);
1568  double v = min;
1569  if (last) n++;
1570  for (int i = 0; i < n; i++) {
1571  array.push_back(v);
1572  v *= binWidth;
1573  }
1574 }
1575 
1585 std::vector<double> AliAnalysisTaskEmcalLight::GenerateLogFixedBinArray(int n, double min, double max, bool last)
1586 {
1587  std::vector<double> array;
1588  GenerateLogFixedBinArray(n, min, max, array, last);
1589  return array;
1590 }
1591 
1592 
1593 TH1* AliAnalysisTaskEmcalLight::GetGeneralTH1(const char* name, bool warn)
1594 {
1595  auto search = fHistograms.find(name);
1596  if (search != fHistograms.end()) {
1597  return search->second;
1598  }
1599  else {
1600  if (warn) AliErrorStream() << "Could not find histogram '" << name << "'" << std::endl;
1601  return nullptr;
1602  }
1603 }
1604 
1605 TH2* AliAnalysisTaskEmcalLight::GetGeneralTH2(const char* name, bool warn)
1606 {
1607  return static_cast<TH2*>(GetGeneralTH1(name, warn));
1608 }
1609 
1610 TProfile* AliAnalysisTaskEmcalLight::GetGeneralTProfile(const char* name, bool warn)
1611 {
1612  return static_cast<TProfile*>(GetGeneralTH1(name, warn));
1613 }
1614 
1616 {
1617  fIsPythia = i;
1618  if (fIsPythia) {
1619  fIsMonteCarlo = kTRUE;
1620  fMCEventHeaderName = "AliGenPythiaEventHeader";
1621  }
1622  else {
1623  if (fMCEventHeaderName == "AliGenPythiaEventHeader") {
1624  fMCEventHeaderName = "";
1625  }
1626  }
1627 }
1628 
1630 {
1631  TClass gen_header_class(name);
1632  if (gen_header_class.InheritsFrom("AliGenEventHeader")) {
1633  fMCEventHeaderName = name;
1634  }
1635  else {
1636  AliWarningStream() << "Class name '" << name << "' does not inherit from 'AliGenEventHeader'. Not setting it." << std::endl;
1637  }
1638 }
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
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
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
Double_t fVertex[3]
!event vertex
Double_t fMaximumEventWeight
Minimum event weight for the related bookkeping histogram.
ECentralityEstimation_t fCentralityEstimation
Centrality estimation.
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.
Bool_t PythiaInfoFromFile(const char *currFile, Float_t &fXsec, Float_t &fTrials, Int_t &pthard)
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)
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 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)