AliPhysics  e6c8d43 (e6c8d43)
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 
522  TProfile* hXsection = GetGeneralTProfile("fHistXsectionNoSel", true);
523  hXsection->SetBinEntries(fPtHardBin + 1, hXsection->GetBinEntries(fPtHardBin + 1) + 1);
524  hXsection->SetBinContent(fPtHardBin + 1, fXsection * hXsection->GetBinEntries(fPtHardBin + 1));
525  }
526 
527  GetGeneralTH1("fHistZVertexNoSel", true)->Fill(fVertex[2]);
528 
529  TH1* hCent = GetGeneralTH1("fHistCentralityNoSel");
530  if (hCent) hCent->Fill(fCent);
531 
532  TH1* hEventPlane = GetGeneralTH1("fHistEventPlaneNoSel");
533  if (hEventPlane) hEventPlane->Fill(fEPV0);
534 
535  TH1* hTriggerClasses = GetGeneralTH1("fHistTriggerClassesNoSel", true);
536  for (auto fired_trg : fFiredTriggerClasses) hTriggerClasses->Fill(fired_trg.c_str(), 1);
537  }
538 
539  return kTRUE;
540 }
541 
562 {
563  if (fInhibit) {
564  AliWarningStream() << "The execution of this task is inhibited. Returning." << std::endl;
565  return;
566  }
567 
568  if (!fLocalInitialized) ExecOnce();
569 
570  if (!fLocalInitialized) return;
571 
572  if (!RetrieveEventObjects()) return;
573 
574  Bool_t eventSelected = IsEventSelected();
575 
577  if (eventSelected) {
578  GetGeneralTH1("fHistEventCount", true)->Fill("Accepted",1);
579  }
580  else {
581  GetGeneralTH1("fHistEventCount", true)->Fill("Rejected",1);
582  }
583 
584  FillGeneralHistograms(kFALSE);
585  if (eventSelected) FillGeneralHistograms(kTRUE);
586  }
587 
588  Bool_t runOk = kFALSE;
589  if (eventSelected || fEventSelectionAfterRun) runOk = Run();
590 
591  if (fCreateHisto && eventSelected && runOk) FillHistograms();
592 
593  if (fCreateHisto && fOutput) {
594  // information for this iteration of the UserExec in the container
595  PostData(1, fOutput);
596  }
597 }
598 
610 Bool_t AliAnalysisTaskEmcalLight::PythiaInfoFromFile(const char* currFile, Float_t &xsec, Float_t &trials, Int_t &pthard)
611 {
612 
613  TString file(currFile);
614  xsec = 0;
615  trials = 1;
616 
617  if (file.Contains(".zip#")) {
618  Ssiz_t pos1 = file.Index("root_archive",12,0,TString::kExact);
619  Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
620  Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
621  file.Replace(pos+1,pos2-pos1,"");
622  } else {
623  // not an archive take the basename....
624  file.ReplaceAll(gSystem->BaseName(file.Data()),"");
625  }
626  AliDebug(1,Form("File name: %s",file.Data()));
627 
628  // Get the pt hard bin
629  TString strPthard(file);
630 
631  strPthard.Remove(strPthard.Last('/'));
632  strPthard.Remove(strPthard.Last('/'));
633  if (strPthard.Contains("AOD")) strPthard.Remove(strPthard.Last('/'));
634  strPthard.Remove(0,strPthard.Last('/')+1);
635  if (strPthard.IsDec()) {
636  pthard = strPthard.Atoi();
637  }
638  else {
639  AliWarning(Form("Could not extract file number from path %s", strPthard.Data()));
640  pthard = -1;
641  }
642 
643  // 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
644  TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root"));
645 
646  if (!fxsec) {
647  // next trial fetch the histgram file
648  fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
649  if (!fxsec) {
650  // not a severe condition but inciate that we have no information
651  return kFALSE;
652  } else {
653  // find the tlist we want to be independtent of the name so use the Tkey
654  TKey* key = static_cast<TKey*>(fxsec->GetListOfKeys()->At(0));
655  if (!key) {
656  fxsec->Close();
657  return kFALSE;
658  }
659  TList *list = dynamic_cast<TList*>(key->ReadObj());
660  if (!list) {
661  fxsec->Close();
662  return kFALSE;
663  }
664  xsec = static_cast<TProfile*>(list->FindObject("h1Xsec"))->GetBinContent(1);
665  trials = static_cast<TH1F*>(list->FindObject("h1Trials"))->GetBinContent(1);
666  fxsec->Close();
667  }
668  } else { // no tree pyxsec.root
669  TTree *xtree = static_cast<TTree*>(fxsec->Get("Xsection"));
670  if (!xtree) {
671  fxsec->Close();
672  return kFALSE;
673  }
674  UInt_t ntrials = 0;
675  Double_t xsection = 0;
676  xtree->SetBranchAddress("xsection",&xsection);
677  xtree->SetBranchAddress("ntrials",&ntrials);
678  xtree->GetEntry(0);
679  trials = ntrials;
680  xsec = xsection;
681  fxsec->Close();
682  }
683  return kTRUE;
684 }
685 
700 {
701  if (!fIsPythia) return kTRUE;
702 
703  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
704  if (!tree) {
705  AliError(Form("%s - UserNotify: No current tree!",GetName()));
706  return kFALSE;
707  }
708 
709  Float_t xsection = 0;
710  Float_t trials = 0;
711  Int_t pthardbin = 0;
712 
713  TFile *curfile = tree->GetCurrentFile();
714  if (!curfile) {
715  AliError(Form("%s - UserNotify: No current file!",GetName()));
716  return kFALSE;
717  }
718 
719  TChain *chain = dynamic_cast<TChain*>(tree);
720  if (chain) tree = chain->GetTree();
721 
722  Int_t nevents = tree->GetEntriesFast();
723 
724  Bool_t res = PythiaInfoFromFile(curfile->GetName(), xsection, trials, pthardbin);
725 
726  fPtHardBin = pthardbin >= 0 ? pthardbin : 0;
727 
728  if (!res) return kTRUE;
729 
731  GetGeneralTH1("fHistTrialsExternalFile", true)->Fill(fPtHardBin, trials);
732  GetGeneralTProfile("fHistXsectionExternalFile", true)->Fill(fPtHardBin, xsection);
733  GetGeneralTH1("fHistEventsExternalFile", true)->Fill(fPtHardBin, nevents);
734  }
735 
736  return kTRUE;
737 }
738 
750 {
751  if (!InputEvent()) {
752  AliError(Form("%s: Could not retrieve event! Returning!", GetName()));
753  return;
754  }
755 
756  if (fNeedEmcalGeom && !fGeom) {
757  fGeom = AliEMCALGeometry::GetInstanceFromRunNumber(InputEvent()->GetRunNumber());
758  if (!fGeom) {
759  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()));
760  return;
761  }
762  }
763 
764  if (fSwitchOffLHC15oFaultyBranches && dynamic_cast<AliAODEvent*>(InputEvent())) {
765  TTree *aodTree = AliAnalysisManager::GetAnalysisManager()->GetTree();
766  aodTree->SetBranchStatus("D0toKpi.fPx", 0);
767  aodTree->SetBranchStatus("D0toKpi.fPy", 0);
768  aodTree->SetBranchStatus("D0toKpi.fPz", 0);
769  aodTree->SetBranchStatus("D0toKpi.fd0", 0);
770  aodTree->SetBranchStatus("Charm3Prong.fPx", 0);
771  aodTree->SetBranchStatus("Charm3Prong.fPy", 0);
772  aodTree->SetBranchStatus("Charm3Prong.fPz", 0);
773  aodTree->SetBranchStatus("Charm3Prong.fd0", 0);
774  aodTree->SetBranchStatus("Dstar.fPx", 0);
775  aodTree->SetBranchStatus("Dstar.fPy", 0);
776  aodTree->SetBranchStatus("Dstar.fPz", 0);
777  aodTree->SetBranchStatus("Dstar.fd0", 0);
778  }
779 
780  //Load all requested track branches - each container knows name already
781  for (auto cont_it : fParticleCollArray) cont_it.second->SetArray(InputEvent());
782 
783  //Load all requested cluster branches - each container knows name already
784  for (auto cont_it : fClusterCollArray) cont_it.second->SetArray(InputEvent());
785 
786  if (!fCaloCellsName.IsNull() && !fCaloCells) {
787  fCaloCells = dynamic_cast<AliVCaloCells*>(InputEvent()->FindListObject(fCaloCellsName));
788  if (!fCaloCells) {
789  AliError(Form("%s: Could not retrieve cells %s!", GetName(), fCaloCellsName.Data()));
790  return;
791  }
792  }
793 
794  if (!fCaloTriggersName.IsNull() && !fCaloTriggers) {
795  fCaloTriggers = dynamic_cast<AliVCaloTrigger*>(InputEvent()->FindListObject(fCaloTriggersName));
796  if (!fCaloTriggers) {
797  AliError(Form("%s: Could not retrieve calo triggers %s!", GetName(), fCaloTriggersName.Data()));
798  return;
799  }
800  }
801 
802  if (!fCaloTriggerPatchInfoName.IsNull() && !fTriggerPatchInfo) {
803  fTriggerPatchInfo = GetArrayFromEvent(fCaloTriggerPatchInfoName.Data(),"AliEMCALTriggerPatchInfo");
804  if (!fTriggerPatchInfo) {
805  AliError(Form("%s: Could not retrieve calo trigger patch info %s!", GetName(), fCaloTriggerPatchInfoName.Data()));
806  return;
807  }
808 
809  }
810 
811  fLocalInitialized = kTRUE;
812 }
813 
820 {
821  if (fForceBeamType != kNA)
822  return fForceBeamType;
823 
824  AliESDEvent *esd = dynamic_cast<AliESDEvent*>(InputEvent());
825  if (esd) {
826  const AliESDRun *run = esd->GetESDRun();
827  TString beamType = run->GetBeamType();
828  if (beamType == "p-p")
829  return kpp;
830  else if (beamType == "A-A")
831  return kAA;
832  else if (beamType == "p-A")
833  return kpA;
834  else
835  return kNA;
836  } else {
837  Int_t runNumber = InputEvent()->GetRunNumber();
838  // All run number ranges taken from the RCT
839  if ((runNumber >= 136833 && runNumber <= 139517) || // LHC10h
840  (runNumber >= 167693 && runNumber <= 170593) || // LHC11h
841  (runNumber >= 244824 && runNumber <= 246994)) { // LHC15o
842  return kAA;
843  } else if ((runNumber >= 188356 && runNumber <= 188366) || // LHC12g
844  (runNumber >= 195164 && runNumber <= 197388) || // LHC13b-f
845  (runNumber >= 265015 && runNumber <= 267166)) { // LHC16q-t
846  return kpA;
847  } else {
848  return kpp;
849  }
850  }
851 }
852 
875 {
876  TH1* hEventRejection = GetGeneralTH1("fHistEventRejection", true);
877 
879  if (fGeneralHistograms) hEventRejection->Fill("PhysSel",1);
880  return kFALSE;
881  }
882 
883  if (!fSelectGeneratorName.IsNull() && !fGeneratorName.IsNull()) {
884  if (!fGeneratorName.Contains(fSelectGeneratorName)) {
885  if (fGeneralHistograms) hEventRejection->Fill("Evt Gen Name",1);
886  return kFALSE;
887  }
888  }
889 
890  Bool_t acceptedTrgClassFound = kFALSE;
891  if (fAcceptedTriggerClasses.size() > 0) {
892  for (auto acc_trg : fAcceptedTriggerClasses) {
893  std::string teststring(acc_trg);
894  bool fullmatch(false);
895  auto posexact = acc_trg.find("EXACT");
896  if(posexact != std::string::npos) {
897  fullmatch = true;
898  teststring.erase(posexact, 5);
899  }
900  for (auto fired_trg : fFiredTriggerClasses) {
901  bool classmatch = fullmatch ? teststring == fired_trg : fired_trg.find(teststring) != std::string::npos;
902  if (classmatch) {
903  acceptedTrgClassFound = kTRUE;
904  break;
905  }
906  }
907  if (acceptedTrgClassFound) break;
908  }
909 
910  if (!acceptedTrgClassFound) {
911  if (fGeneralHistograms) hEventRejection->Fill("Trg class (acc)",1);
912  return kFALSE;
913  }
914  }
915 
916  if (fRejectedTriggerClasses.size() > 0) {
917  for (auto rej_trg : fRejectedTriggerClasses) {
918  std::string teststring(rej_trg);
919  bool fullmatch(false);
920  auto posexact = rej_trg.find("EXACT");
921  if(posexact != std::string::npos) {
922  fullmatch = true;
923  teststring.erase(posexact, 5);
924  }
925  for (auto fired_trg : fFiredTriggerClasses) {
926  bool classmatch = fullmatch ? teststring == fired_trg : fired_trg.find(teststring) != std::string::npos;
927  if (classmatch) {
928  if (fGeneralHistograms) hEventRejection->Fill("Trg class (rej)",1);
929  return kFALSE;
930  }
931  }
932  }
933  }
934 
935  if (fMinCent < fMaxCent && fMaxCent > 0) {
936  if (fCent < fMinCent || fCent > fMaxCent) {
937  if (fGeneralHistograms) hEventRejection->Fill("Cent",1);
938  return kFALSE;
939  }
940  }
941 
942  if (fNVertCont < fMinNVertCont) {
943  if (fGeneralHistograms) hEventRejection->Fill("vertex contr.",1);
944  return kFALSE;
945  }
946 
947  if (fMinVz < fMaxVz) {
948  if (fVertex[2] < fMinVz || fVertex[2] > fMaxVz) {
949  if (fGeneralHistograms) hEventRejection->Fill("Vz",1);
950  return kFALSE;
951  }
952  }
953 
954  if (fMaxVzDiff >= 0) {
955  if (fNVertSPDCont > 0) {
956  Double_t vzSPD = fVertexSPD[2];
957  Double_t dvertex = TMath::Abs(fVertex[2] - vzSPD);
958  //if difference larger than fZvertexDiff
959  if (dvertex > fMaxVzDiff) {
960  if (fGeneralHistograms) hEventRejection->Fill("VzSPD",1);
961  return kFALSE;
962  }
963  }
964  }
965 
966  if (fMinPtHard >= 0 && fPtHard < fMinPtHard) {
967  if (fGeneralHistograms) hEventRejection->Fill("SelPtHardBin",1);
968  return kFALSE;
969  }
970 
971  if (fMaxPtHard >= 0 && fPtHard >= fMaxPtHard) {
972  if (fGeneralHistograms) hEventRejection->Fill("SelPtHardBin",1);
973  return kFALSE;
974  }
975 
977  if (fGeneralHistograms) hEventRejection->Fill("SelPtHardBin",1);
978  return kFALSE;
979  }
980 
981  // Reject filter for MC data
982  if (!CheckMCOutliers()) {
983  if (fGeneralHistograms) hEventRejection->Fill("MCOutlier",1);
984  return kFALSE;
985  }
986 
987  return kTRUE;
988 }
989 
998 TClonesArray *AliAnalysisTaskEmcalLight::GetArrayFromEvent(const char *name, const char *clname)
999 {
1000  TClonesArray *arr = 0;
1001  TString sname(name);
1002  if (!sname.IsNull()) {
1003  arr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(sname));
1004  if (!arr) {
1005  AliWarning(Form("%s: Could not retrieve array with name %s!", GetName(), name));
1006  return 0;
1007  }
1008  } else {
1009  return 0;
1010  }
1011 
1012  if (!clname)
1013  return arr;
1014 
1015  TString objname(arr->GetClass()->GetName());
1016  TClass cls(objname);
1017  if (!cls.InheritsFrom(clname)) {
1018  AliWarning(Form("%s: Objects of type %s in %s are not inherited from %s!",
1019  GetName(), cls.GetName(), name, clname));
1020  return 0;
1021  }
1022  return arr;
1023 }
1024 
1030 {
1031  fVertex[0] = 0;
1032  fVertex[1] = 0;
1033  fVertex[2] = 0;
1034  fNVertCont = 0;
1035 
1036  fVertexSPD[0] = 0;
1037  fVertexSPD[1] = 0;
1038  fVertexSPD[2] = 0;
1039  fNVertSPDCont = 0;
1040 
1041  fFiredTriggerClasses.clear();
1042  std::stringstream firedClasses(InputEvent()->GetFiredTriggerClasses().Data());
1043  while (firedClasses.good()) {
1044  std::string trgClass;
1045  firedClasses >> trgClass;
1046  if (!trgClass.empty()) fFiredTriggerClasses.push_back(trgClass);
1047  }
1048 
1049  if (fDataType == kESD) {
1050  fFiredTriggerBitMap = static_cast<AliInputEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler())->IsEventSelected();
1051  }
1052  else {
1053  fFiredTriggerBitMap = static_cast<AliVAODHeader*>(InputEvent()->GetHeader())->GetOfflineTrigger();
1054  }
1055 
1056  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
1057  if (vert) {
1058  vert->GetXYZ(fVertex);
1059  fNVertCont = vert->GetNContributors();
1060  }
1061 
1062  const AliVVertex *vertSPD = InputEvent()->GetPrimaryVertexSPD();
1063  if (vertSPD) {
1064  vertSPD->GetXYZ(fVertexSPD);
1065  fNVertSPDCont = vertSPD->GetNContributors();
1066  }
1067 
1068  fBeamType = GetBeamType();
1069 
1070  fCent = 99;
1071  fCentBin = -1;
1072  fEPV0 = -999;
1073  fEPV0A = -999;
1074  fEPV0C = -999;
1075 
1077  // New centrality estimation (AliMultSelection)
1078  // See https://twiki.cern.ch/twiki/bin/viewauth/ALICE/AliMultSelectionCalibStatus for calibration status period-by-period)
1079  AliMultSelection *MultSelection = static_cast<AliMultSelection*>(InputEvent()->FindListObject("MultSelection"));
1080  if (MultSelection) {
1081  fCent = MultSelection->GetMultiplicityPercentile(fCentEst.Data());
1082  }
1083  else {
1084  AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
1085  }
1086  }
1087  else if (fCentralityEstimation == kOldCentrality) {
1088  // Old centrality estimation (AliCentrality, works only on Run-1 PbPb and pPb)
1089  AliCentrality *aliCent = InputEvent()->GetCentrality();
1090  if (aliCent) {
1091  fCent = aliCent->GetCentralityPercentile(fCentEst.Data());
1092  }
1093  else {
1094  AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
1095  }
1096  }
1097  if (!fCentBins.empty() && fCentralityEstimation != kNoCentrality) {
1098  for (auto cent_it = fCentBins.begin(); cent_it != fCentBins.end() - 1; cent_it++) {
1099  if (fCent >= *cent_it && fCent < *(cent_it+1)) fCentBin = cent_it - fCentBins.begin();
1100  }
1101  }
1102  else {
1103  fCentBin = 0;
1104  }
1105 
1106  if (fBeamType == kAA || fBeamType == kpA ) {
1107  AliEventplane *aliEP = InputEvent()->GetEventplane();
1108  if (aliEP) {
1109  fEPV0 = aliEP->GetEventplane("V0" ,InputEvent());
1110  fEPV0A = aliEP->GetEventplane("V0A",InputEvent());
1111  fEPV0C = aliEP->GetEventplane("V0C",InputEvent());
1112  } else {
1113  AliWarning(Form("%s: Could not retrieve event plane information!", GetName()));
1114  }
1115  }
1116 
1117  if (fIsMonteCarlo && MCEvent()) {
1118  AliGenEventHeader* header = MCEvent()->GenEventHeader();
1119  if (fMCEventHeaderName.IsNull()) {
1120  fMCHeader = header;
1121  }
1122  else {
1123  if (header->InheritsFrom(fMCEventHeaderName)) {
1124  fMCHeader = header;
1125  }
1126  else if (header->InheritsFrom("AliGenCocktailEventHeader")) {
1127  AliGenCocktailEventHeader* cocktailHeader = static_cast<AliGenCocktailEventHeader*>(header);
1128  TList* headers = cocktailHeader->GetHeaders();
1129  for (auto obj : *headers) { // @suppress("Symbol is not resolved")
1130  if (obj->InheritsFrom(fMCEventHeaderName)){
1131  fMCHeader = static_cast<AliGenEventHeader*>(obj);
1132  break;
1133  }
1134  }
1135  }
1136  }
1137  if (fMCHeader) {
1138  fEventWeight = fMCHeader->EventWeight();
1139  if (fIsPythia) {
1140  fPythiaHeader = static_cast<AliGenPythiaEventHeader*>(fMCHeader);
1141  fPtHard = fPythiaHeader->GetPtHard();
1142  fXsection = fPythiaHeader->GetXsection();
1143  fNTrials = fPythiaHeader->Trials();
1144  }
1145  }
1146  }
1147 
1148  for (auto cont_it : fParticleCollArray) cont_it.second->NextEvent(InputEvent());
1149  for (auto cont_it : fClusterCollArray) cont_it.second->NextEvent(InputEvent());
1150 
1151  return kTRUE;
1152 }
1153 
1162 AliParticleContainer* AliAnalysisTaskEmcalLight::AddParticleContainer(std::string branchName, std::string contName)
1163 {
1164  if (branchName.size() == 0) return 0;
1165 
1166  AliParticleContainer* cont = 0;
1167 
1168  if (branchName == "tracks" || branchName == "Tracks") cont = new AliTrackContainer(branchName.c_str());
1169  else if (branchName == "mcparticles") cont = new AliMCParticleContainer(branchName.c_str());
1170  else cont = new AliParticleContainer(branchName.c_str());
1171 
1172  if (contName.size() > 0) cont->SetName(contName.c_str());
1173 
1174  AdoptParticleContainer(cont);
1175 
1176  return cont;
1177 }
1178 
1187 AliClusterContainer* AliAnalysisTaskEmcalLight::AddClusterContainer(std::string branchName, std::string contName)
1188 {
1189  if (branchName.size() == 0) return 0;
1190 
1191  AliClusterContainer* cont = new AliClusterContainer(branchName.c_str());
1192 
1193  if (contName.size() > 0) cont->SetName(contName.c_str());
1194 
1195  AdoptClusterContainer(cont);
1196 
1197  return cont;
1198 }
1199 
1206 {
1207  std::map<std::string, AliParticleContainer*>::const_iterator cont_it = fParticleCollArray.find(name);
1208  if (cont_it != fParticleCollArray.end()) return cont_it->second;
1209  else return nullptr;
1210 }
1211 
1218 {
1219  std::map<std::string, AliClusterContainer*>::const_iterator cont_it = fClusterCollArray.find(name);
1220  if (cont_it != fClusterCollArray.end()) return cont_it->second;
1221  else return nullptr;
1222 }
1223 
1230 {
1231  if (!(InputEvent()->FindListObject(obj->GetName()))) {
1232  InputEvent()->AddObject(obj);
1233  }
1234  else {
1235  if (!attempt) {
1236  AliFatal(Form("%s: Container with name %s already present. Aborting", GetName(), obj->GetName()));
1237  }
1238  }
1239 }
1240 
1249 {
1250 
1251  if (!fGeom) {
1252  AliWarning(Form("%s - AliAnalysisTaskEmcalBase::IsTrackInEmcalAcceptance - Geometry is not available!", GetName()));
1253  return kFALSE;
1254  }
1255 
1256  Double_t minPhi = fGeom->GetArm1PhiMin() - edges;
1257  Double_t maxPhi = fGeom->GetArm1PhiMax() + edges;
1258 
1259  if (part->Phi() > minPhi && part->Phi() < maxPhi) {
1260  return kTRUE;
1261  }
1262  else {
1263  return kFALSE;
1264  }
1265 }
1266 
1268 {
1269  axis->SetBinLabel(1, "NullObject");
1270  axis->SetBinLabel(2, "Pt");
1271  axis->SetBinLabel(3, "Acceptance");
1272  axis->SetBinLabel(4, "MCLabel");
1273  axis->SetBinLabel(5, "BitMap");
1274  axis->SetBinLabel(6, "HF cut");
1275  axis->SetBinLabel(7, "Bit6");
1276  axis->SetBinLabel(8, "NotHybridTrack");
1277  axis->SetBinLabel(9, "MCFlag");
1278  axis->SetBinLabel(10, "MCGenerator");
1279  axis->SetBinLabel(11, "ChargeCut");
1280  axis->SetBinLabel(12, "MinDistanceTPCSectorEdge");
1281  axis->SetBinLabel(13, "Bit12");
1282  axis->SetBinLabel(14, "IsEMCal");
1283  axis->SetBinLabel(15, "Time");
1284  axis->SetBinLabel(16, "Energy");
1285  axis->SetBinLabel(17, "ExoticCut");
1286  axis->SetBinLabel(18, "Bit17");
1287  axis->SetBinLabel(19, "Area");
1288  axis->SetBinLabel(20, "AreaEmc");
1289  axis->SetBinLabel(21, "ZLeadingCh");
1290  axis->SetBinLabel(22, "ZLeadingEmc");
1291  axis->SetBinLabel(23, "NEF");
1292  axis->SetBinLabel(24, "MinLeadPt");
1293  axis->SetBinLabel(25, "MaxTrackPt");
1294  axis->SetBinLabel(26, "MaxClusterPt");
1295  axis->SetBinLabel(27, "Flavour");
1296  axis->SetBinLabel(28, "TagStatus");
1297  axis->SetBinLabel(29, "MinNConstituents");
1298  axis->SetBinLabel(30, "Bit29");
1299  axis->SetBinLabel(31, "Bit30");
1300  axis->SetBinLabel(32, "Bit31");
1301 }
1302 
1309 Double_t AliAnalysisTaskEmcalLight::GetParallelFraction(AliVParticle* part1, AliVParticle* part2)
1310 {
1311  TVector3 vect1(part1->Px(), part1->Py(), part1->Pz());
1312  TVector3 vect2(part2->Px(), part2->Py(), part2->Pz());
1313  Double_t z = (vect1 * vect2) / (vect2 * vect2);
1314  return z;
1315 }
1316 
1323 Double_t AliAnalysisTaskEmcalLight::GetParallelFraction(const TVector3& vect1, AliVParticle* part2)
1324 {
1325  TVector3 vect2(part2->Px(), part2->Py(), part2->Pz());
1326  Double_t z = (vect1 * vect2) / (vect2 * vect2);
1327  return z;
1328 }
1329 
1338 void AliAnalysisTaskEmcalLight::GetEtaPhiDiff(const AliVTrack *t, const AliVCluster *v, Double_t &phidiff, Double_t &etadiff)
1339 {
1340  phidiff = 999;
1341  etadiff = 999;
1342 
1343  if (!t||!v) return;
1344 
1345  Double_t veta = t->GetTrackEtaOnEMCal();
1346  Double_t vphi = t->GetTrackPhiOnEMCal();
1347 
1348  Float_t pos[3] = {0};
1349  v->GetPosition(pos);
1350  TVector3 cpos(pos);
1351  Double_t ceta = cpos.Eta();
1352  Double_t cphi = cpos.Phi();
1353  etadiff=veta-ceta;
1354  phidiff=TVector2::Phi_mpi_pi(vphi-cphi);
1355 }
1356 
1363 {
1364  Byte_t ret = 0;
1365  if (t->TestBit(BIT(22)) && !t->TestBit(BIT(23)))
1366  ret = 1;
1367  else if (!t->TestBit(BIT(22)) && t->TestBit(BIT(23)))
1368  ret = 2;
1369  else if (t->TestBit(BIT(22)) && t->TestBit(BIT(23)))
1370  ret = 3;
1371  return ret;
1372 }
1373 
1383 Byte_t AliAnalysisTaskEmcalLight::GetTrackType(const AliAODTrack *aodTrack, UInt_t filterBit1, UInt_t filterBit2)
1384 {
1385 
1386  Int_t res = 0;
1387 
1388  if (aodTrack->TestFilterBit(filterBit1)) {
1389  res = 0;
1390  }
1391  else if (aodTrack->TestFilterBit(filterBit2)) {
1392  if ((aodTrack->GetStatus()&AliVTrack::kITSrefit)!=0) {
1393  res = 1;
1394  }
1395  else {
1396  res = 2;
1397  }
1398  }
1399  else {
1400  res = 3;
1401  }
1402 
1403  return res;
1404 }
1405 
1412 {
1413  EBeamType_t b = kpp;
1414  if ((runnumber >= 136833 && runnumber <= 139517) || // LHC10h Run-1 (Pb-Pb)
1415  (runnumber >= 167693 && runnumber <= 170593) || // LHC11h Run-1 (Pb-Pb)
1416  (runnumber >= 244824 && runnumber <= 246994)) { // LHC15o Run-2 (Pb-Pb)
1417  b = kAA;
1418  }
1419  else if ((runnumber > 188356 && runnumber <= 188503) || // LHC12g Run-1 (p-Pb pilot)
1420  (runnumber >= 195164 && runnumber <= 197388) || // LHC13b,c,d,e,f Run-1 (p-Pb)
1421  (runnumber >= 265077 && runnumber <= 267166)) { // LHC16 Run-2 (p-Pb)
1422  b = kpA;
1423  }
1424  return b;
1425 }
1426 
1433 {
1434  if (!fPythiaHeader || !fMCRejectFilter) return kTRUE;
1435 
1436  // Condition 1: Pythia jet / pT-hard > factor
1437  if (fPtHardAndJetPtFactor > 0.) {
1438  AliTLorentzVector jet;
1439 
1440  Int_t nTriggerJets = fPythiaHeader->NTriggerJets();
1441 
1442  AliDebug(1,Form("Njets: %d, pT Hard %f",nTriggerJets, fPtHard));
1443 
1444  Float_t tmpjet[]={0,0,0,0};
1445  for (Int_t ijet = 0; ijet< nTriggerJets; ijet++) {
1446  fPythiaHeader->TriggerJet(ijet, tmpjet);
1447 
1448  jet.SetPxPyPzE(tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3]);
1449 
1450  AliDebug(1,Form("jet %d; pycell jet pT %f",ijet, jet.Pt()));
1451 
1452  //Compare jet pT and pt Hard
1453  if (jet.Pt() > fPtHardAndJetPtFactor * fPtHard) {
1454  AliInfo(Form("Reject jet event with : pT Hard %2.2f, pycell jet pT %2.2f, rejection factor %1.1f\n", fPtHard, jet.Pt(), fPtHardAndJetPtFactor));
1455  return kFALSE;
1456  }
1457  }
1458  }
1459  // end condition 1
1460 
1461  // Condition 2 : Reconstructed EMCal cluster pT / pT-hard > factor
1462  if (fPtHardAndClusterPtFactor > 0.) {
1463  AliClusterContainer* mccluscont = GetClusterContainer(0);
1464  if ((Bool_t)mccluscont) {
1465  for (auto cluster : mccluscont->all()) {// Not cuts applied ; use accept for cuts
1466  Float_t ecluster = cluster->E();
1467 
1468  if (ecluster > (fPtHardAndClusterPtFactor * fPtHard)) {
1469  AliInfo(Form("Reject : ecluster %2.2f, calo %d, factor %2.2f, ptHard %f",ecluster,cluster->GetType(),fPtHardAndClusterPtFactor,fPtHard));
1470  return kFALSE;
1471  }
1472  }
1473  }
1474  }
1475  // end condition 2
1476 
1477  // condition 3 : Reconstructed track pT / pT-hard >factor
1478  if (fPtHardAndTrackPtFactor > 0.) {
1479  AliMCParticleContainer* mcpartcont = dynamic_cast<AliMCParticleContainer*>(GetParticleContainer(0));
1480  if ((Bool_t)mcpartcont) {
1481  for (auto mctrack : mcpartcont->all()) {// Not cuts applied ; use accept for cuts
1482  Float_t trackpt = mctrack->Pt();
1483  if (trackpt > (fPtHardAndTrackPtFactor * fPtHard) ) {
1484  AliInfo(Form("Reject : track %2.2f, factor %2.2f, ptHard %f", trackpt, fPtHardAndTrackPtFactor, fPtHard));
1485  return kFALSE;
1486  }
1487  }
1488  }
1489  }
1490  // end condition 3
1491 
1492  return kTRUE;
1493 }
1494 
1504 {
1505  Double_t dphi = -999;
1506  const Double_t tpi = TMath::TwoPi();
1507 
1508  if (phia < 0) phia += tpi;
1509  else if (phia > tpi) phia -= tpi;
1510  if (phib < 0) phib += tpi;
1511  else if (phib > tpi) phib -= tpi;
1512  dphi = phib - phia;
1513  if (dphi < rangeMin) dphi += tpi;
1514  else if (dphi > rangeMax) dphi -= tpi;
1515 
1516  return dphi;
1517 }
1518 
1528 void AliAnalysisTaskEmcalLight::GenerateFixedBinArray(int n, double min, double max, std::vector<double>& array, bool last)
1529 {
1530  double binWidth = (max - min) / n;
1531  double v = min;
1532  if (last) n++;
1533  for (int i = 0; i < n; i++) {
1534  array.push_back(v);
1535  v += binWidth;
1536  }
1537 }
1538 
1548 std::vector<double> AliAnalysisTaskEmcalLight::GenerateFixedBinArray(int n, double min, double max, bool last)
1549 {
1550  std::vector<double> array;
1551  GenerateFixedBinArray(n, min, max, array, last);
1552  return array;
1553 }
1554 
1564 void AliAnalysisTaskEmcalLight::GenerateLogFixedBinArray(int n, double min, double max, std::vector<double>& array, bool last)
1565 {
1566  if (min <= 0 || max < min) {
1567  AliErrorClassStream() << "Cannot generate a log scale fixed-bin array with limits " << min << ", " << max << std::endl;
1568  return;
1569  }
1570  double binWidth = std::pow(max / min, 1.0 / n);
1571  double v = min;
1572  if (last) n++;
1573  for (int i = 0; i < n; i++) {
1574  array.push_back(v);
1575  v *= binWidth;
1576  }
1577 }
1578 
1588 std::vector<double> AliAnalysisTaskEmcalLight::GenerateLogFixedBinArray(int n, double min, double max, bool last)
1589 {
1590  std::vector<double> array;
1591  GenerateLogFixedBinArray(n, min, max, array, last);
1592  return array;
1593 }
1594 
1595 
1596 TH1* AliAnalysisTaskEmcalLight::GetGeneralTH1(const char* name, bool warn)
1597 {
1598  auto search = fHistograms.find(name);
1599  if (search != fHistograms.end()) {
1600  return search->second;
1601  }
1602  else {
1603  if (warn) AliErrorStream() << "Could not find histogram '" << name << "'" << std::endl;
1604  return nullptr;
1605  }
1606 }
1607 
1608 TH2* AliAnalysisTaskEmcalLight::GetGeneralTH2(const char* name, bool warn)
1609 {
1610  return static_cast<TH2*>(GetGeneralTH1(name, warn));
1611 }
1612 
1613 TProfile* AliAnalysisTaskEmcalLight::GetGeneralTProfile(const char* name, bool warn)
1614 {
1615  return static_cast<TProfile*>(GetGeneralTH1(name, warn));
1616 }
1617 
1619 {
1620  fIsPythia = i;
1621  if (fIsPythia) {
1622  fIsMonteCarlo = kTRUE;
1623  fMCEventHeaderName = "AliGenPythiaEventHeader";
1624  }
1625  else {
1626  if (fMCEventHeaderName == "AliGenPythiaEventHeader") {
1627  fMCEventHeaderName = "";
1628  }
1629  }
1630 }
1631 
1633 {
1634  TClass gen_header_class(name);
1635  if (gen_header_class.InheritsFrom("AliGenEventHeader")) {
1636  fMCEventHeaderName = name;
1637  }
1638  else {
1639  AliWarningStream() << "Class name '" << name << "' does not inherit from 'AliGenEventHeader'. Not setting it." << std::endl;
1640  }
1641 }
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)