AliPhysics  2853087 (2853087)
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  for (auto fired_trg : fFiredTriggerClasses) {
894  if (fired_trg.find(acc_trg) != std::string::npos) {
895  acceptedTrgClassFound = kTRUE;
896  break;
897  }
898  }
899  if (acceptedTrgClassFound) break;
900  }
901 
902  if (!acceptedTrgClassFound) {
903  if (fGeneralHistograms) hEventRejection->Fill("Trg class (acc)",1);
904  return kFALSE;
905  }
906  }
907 
908  if (fRejectedTriggerClasses.size() > 0) {
909  for (auto rej_trg : fRejectedTriggerClasses) {
910  for (auto fired_trg : fFiredTriggerClasses) {
911  if (fired_trg.find(rej_trg) != std::string::npos) {
912  if (fGeneralHistograms) hEventRejection->Fill("Trg class (rej)",1);
913  return kFALSE;
914  }
915  }
916  }
917  }
918 
919  if (fMinCent < fMaxCent && fMaxCent > 0) {
920  if (fCent < fMinCent || fCent > fMaxCent) {
921  if (fGeneralHistograms) hEventRejection->Fill("Cent",1);
922  return kFALSE;
923  }
924  }
925 
926  if (fNVertCont < fMinNVertCont) {
927  if (fGeneralHistograms) hEventRejection->Fill("vertex contr.",1);
928  return kFALSE;
929  }
930 
931  if (fMinVz < fMaxVz) {
932  if (fVertex[2] < fMinVz || fVertex[2] > fMaxVz) {
933  if (fGeneralHistograms) hEventRejection->Fill("Vz",1);
934  return kFALSE;
935  }
936  }
937 
938  if (fMaxVzDiff >= 0) {
939  if (fNVertSPDCont > 0) {
940  Double_t vzSPD = fVertexSPD[2];
941  Double_t dvertex = TMath::Abs(fVertex[2] - vzSPD);
942  //if difference larger than fZvertexDiff
943  if (dvertex > fMaxVzDiff) {
944  if (fGeneralHistograms) hEventRejection->Fill("VzSPD",1);
945  return kFALSE;
946  }
947  }
948  }
949 
950  if (fMinPtHard >= 0 && fPtHard < fMinPtHard) {
951  if (fGeneralHistograms) hEventRejection->Fill("SelPtHardBin",1);
952  return kFALSE;
953  }
954 
955  if (fMaxPtHard >= 0 && fPtHard >= fMaxPtHard) {
956  if (fGeneralHistograms) hEventRejection->Fill("SelPtHardBin",1);
957  return kFALSE;
958  }
959 
961  if (fGeneralHistograms) hEventRejection->Fill("SelPtHardBin",1);
962  return kFALSE;
963  }
964 
965  // Reject filter for MC data
966  if (!CheckMCOutliers()) {
967  if (fGeneralHistograms) hEventRejection->Fill("MCOutlier",1);
968  return kFALSE;
969  }
970 
971  return kTRUE;
972 }
973 
982 TClonesArray *AliAnalysisTaskEmcalLight::GetArrayFromEvent(const char *name, const char *clname)
983 {
984  TClonesArray *arr = 0;
985  TString sname(name);
986  if (!sname.IsNull()) {
987  arr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(sname));
988  if (!arr) {
989  AliWarning(Form("%s: Could not retrieve array with name %s!", GetName(), name));
990  return 0;
991  }
992  } else {
993  return 0;
994  }
995 
996  if (!clname)
997  return arr;
998 
999  TString objname(arr->GetClass()->GetName());
1000  TClass cls(objname);
1001  if (!cls.InheritsFrom(clname)) {
1002  AliWarning(Form("%s: Objects of type %s in %s are not inherited from %s!",
1003  GetName(), cls.GetName(), name, clname));
1004  return 0;
1005  }
1006  return arr;
1007 }
1008 
1014 {
1015  fVertex[0] = 0;
1016  fVertex[1] = 0;
1017  fVertex[2] = 0;
1018  fNVertCont = 0;
1019 
1020  fVertexSPD[0] = 0;
1021  fVertexSPD[1] = 0;
1022  fVertexSPD[2] = 0;
1023  fNVertSPDCont = 0;
1024 
1025  fFiredTriggerClasses.clear();
1026  std::stringstream firedClasses(InputEvent()->GetFiredTriggerClasses().Data());
1027  while (firedClasses.good()) {
1028  std::string trgClass;
1029  firedClasses >> trgClass;
1030  if (!trgClass.empty()) fFiredTriggerClasses.push_back(trgClass);
1031  }
1032 
1033  if (fDataType == kESD) {
1034  fFiredTriggerBitMap = static_cast<AliInputEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler())->IsEventSelected();
1035  }
1036  else {
1037  fFiredTriggerBitMap = static_cast<AliVAODHeader*>(InputEvent()->GetHeader())->GetOfflineTrigger();
1038  }
1039 
1040  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
1041  if (vert) {
1042  vert->GetXYZ(fVertex);
1043  fNVertCont = vert->GetNContributors();
1044  }
1045 
1046  const AliVVertex *vertSPD = InputEvent()->GetPrimaryVertexSPD();
1047  if (vertSPD) {
1048  vertSPD->GetXYZ(fVertexSPD);
1049  fNVertSPDCont = vertSPD->GetNContributors();
1050  }
1051 
1052  fBeamType = GetBeamType();
1053 
1054  fCent = 99;
1055  fCentBin = -1;
1056  fEPV0 = -999;
1057  fEPV0A = -999;
1058  fEPV0C = -999;
1059 
1061  // New centrality estimation (AliMultSelection)
1062  // See https://twiki.cern.ch/twiki/bin/viewauth/ALICE/AliMultSelectionCalibStatus for calibration status period-by-period)
1063  AliMultSelection *MultSelection = static_cast<AliMultSelection*>(InputEvent()->FindListObject("MultSelection"));
1064  if (MultSelection) {
1065  fCent = MultSelection->GetMultiplicityPercentile(fCentEst.Data());
1066  }
1067  else {
1068  AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
1069  }
1070  }
1071  else if (fCentralityEstimation == kOldCentrality) {
1072  // Old centrality estimation (AliCentrality, works only on Run-1 PbPb and pPb)
1073  AliCentrality *aliCent = InputEvent()->GetCentrality();
1074  if (aliCent) {
1075  fCent = aliCent->GetCentralityPercentile(fCentEst.Data());
1076  }
1077  else {
1078  AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
1079  }
1080  }
1081  if (!fCentBins.empty() && fCentralityEstimation != kNoCentrality) {
1082  for (auto cent_it = fCentBins.begin(); cent_it != fCentBins.end() - 1; cent_it++) {
1083  if (fCent >= *cent_it && fCent < *(cent_it+1)) fCentBin = cent_it - fCentBins.begin();
1084  }
1085  }
1086  else {
1087  fCentBin = 0;
1088  }
1089 
1090  if (fBeamType == kAA || fBeamType == kpA ) {
1091  AliEventplane *aliEP = InputEvent()->GetEventplane();
1092  if (aliEP) {
1093  fEPV0 = aliEP->GetEventplane("V0" ,InputEvent());
1094  fEPV0A = aliEP->GetEventplane("V0A",InputEvent());
1095  fEPV0C = aliEP->GetEventplane("V0C",InputEvent());
1096  } else {
1097  AliWarning(Form("%s: Could not retrieve event plane information!", GetName()));
1098  }
1099  }
1100 
1101  if (fIsMonteCarlo && MCEvent()) {
1102  AliGenEventHeader* header = MCEvent()->GenEventHeader();
1103  if (fMCEventHeaderName.IsNull()) {
1104  fMCHeader = header;
1105  }
1106  else {
1107  if (header->InheritsFrom(fMCEventHeaderName)) {
1108  fMCHeader = header;
1109  }
1110  else if (header->InheritsFrom("AliGenCocktailEventHeader")) {
1111  AliGenCocktailEventHeader* cocktailHeader = static_cast<AliGenCocktailEventHeader*>(header);
1112  TList* headers = cocktailHeader->GetHeaders();
1113  for (auto obj : *headers) { // @suppress("Symbol is not resolved")
1114  if (obj->InheritsFrom(fMCEventHeaderName)){
1115  fMCHeader = static_cast<AliGenEventHeader*>(obj);
1116  break;
1117  }
1118  }
1119  }
1120  }
1121  if (fMCHeader) {
1122  fEventWeight = fMCHeader->EventWeight();
1123  if (fIsPythia) {
1124  fPythiaHeader = static_cast<AliGenPythiaEventHeader*>(fMCHeader);
1125  fPtHard = fPythiaHeader->GetPtHard();
1126  fXsection = fPythiaHeader->GetXsection();
1127  fNTrials = fPythiaHeader->Trials();
1128  }
1129  }
1130  }
1131 
1132  for (auto cont_it : fParticleCollArray) cont_it.second->NextEvent(InputEvent());
1133  for (auto cont_it : fClusterCollArray) cont_it.second->NextEvent(InputEvent());
1134 
1135  return kTRUE;
1136 }
1137 
1146 AliParticleContainer* AliAnalysisTaskEmcalLight::AddParticleContainer(std::string branchName, std::string contName)
1147 {
1148  if (branchName.size() == 0) return 0;
1149 
1150  AliParticleContainer* cont = 0;
1151 
1152  if (branchName == "tracks" || branchName == "Tracks") cont = new AliTrackContainer(branchName.c_str());
1153  else if (branchName == "mcparticles") cont = new AliMCParticleContainer(branchName.c_str());
1154  else cont = new AliParticleContainer(branchName.c_str());
1155 
1156  if (contName.size() > 0) cont->SetName(contName.c_str());
1157 
1158  AdoptParticleContainer(cont);
1159 
1160  return cont;
1161 }
1162 
1171 AliClusterContainer* AliAnalysisTaskEmcalLight::AddClusterContainer(std::string branchName, std::string contName)
1172 {
1173  if (branchName.size() == 0) return 0;
1174 
1175  AliClusterContainer* cont = new AliClusterContainer(branchName.c_str());
1176 
1177  if (contName.size() > 0) cont->SetName(contName.c_str());
1178 
1179  AdoptClusterContainer(cont);
1180 
1181  return cont;
1182 }
1183 
1190 {
1191  std::map<std::string, AliParticleContainer*>::const_iterator cont_it = fParticleCollArray.find(name);
1192  if (cont_it != fParticleCollArray.end()) return cont_it->second;
1193  else return nullptr;
1194 }
1195 
1202 {
1203  std::map<std::string, AliClusterContainer*>::const_iterator cont_it = fClusterCollArray.find(name);
1204  if (cont_it != fClusterCollArray.end()) return cont_it->second;
1205  else return nullptr;
1206 }
1207 
1214 {
1215  if (!(InputEvent()->FindListObject(obj->GetName()))) {
1216  InputEvent()->AddObject(obj);
1217  }
1218  else {
1219  if (!attempt) {
1220  AliFatal(Form("%s: Container with name %s already present. Aborting", GetName(), obj->GetName()));
1221  }
1222  }
1223 }
1224 
1233 {
1234 
1235  if (!fGeom) {
1236  AliWarning(Form("%s - AliAnalysisTaskEmcalBase::IsTrackInEmcalAcceptance - Geometry is not available!", GetName()));
1237  return kFALSE;
1238  }
1239 
1240  Double_t minPhi = fGeom->GetArm1PhiMin() - edges;
1241  Double_t maxPhi = fGeom->GetArm1PhiMax() + edges;
1242 
1243  if (part->Phi() > minPhi && part->Phi() < maxPhi) {
1244  return kTRUE;
1245  }
1246  else {
1247  return kFALSE;
1248  }
1249 }
1250 
1252 {
1253  axis->SetBinLabel(1, "NullObject");
1254  axis->SetBinLabel(2, "Pt");
1255  axis->SetBinLabel(3, "Acceptance");
1256  axis->SetBinLabel(4, "MCLabel");
1257  axis->SetBinLabel(5, "BitMap");
1258  axis->SetBinLabel(6, "HF cut");
1259  axis->SetBinLabel(7, "Bit6");
1260  axis->SetBinLabel(8, "NotHybridTrack");
1261  axis->SetBinLabel(9, "MCFlag");
1262  axis->SetBinLabel(10, "MCGenerator");
1263  axis->SetBinLabel(11, "ChargeCut");
1264  axis->SetBinLabel(12, "MinDistanceTPCSectorEdge");
1265  axis->SetBinLabel(13, "Bit12");
1266  axis->SetBinLabel(14, "IsEMCal");
1267  axis->SetBinLabel(15, "Time");
1268  axis->SetBinLabel(16, "Energy");
1269  axis->SetBinLabel(17, "ExoticCut");
1270  axis->SetBinLabel(18, "Bit17");
1271  axis->SetBinLabel(19, "Area");
1272  axis->SetBinLabel(20, "AreaEmc");
1273  axis->SetBinLabel(21, "ZLeadingCh");
1274  axis->SetBinLabel(22, "ZLeadingEmc");
1275  axis->SetBinLabel(23, "NEF");
1276  axis->SetBinLabel(24, "MinLeadPt");
1277  axis->SetBinLabel(25, "MaxTrackPt");
1278  axis->SetBinLabel(26, "MaxClusterPt");
1279  axis->SetBinLabel(27, "Flavour");
1280  axis->SetBinLabel(28, "TagStatus");
1281  axis->SetBinLabel(29, "MinNConstituents");
1282  axis->SetBinLabel(30, "Bit29");
1283  axis->SetBinLabel(31, "Bit30");
1284  axis->SetBinLabel(32, "Bit31");
1285 }
1286 
1293 Double_t AliAnalysisTaskEmcalLight::GetParallelFraction(AliVParticle* part1, AliVParticle* part2)
1294 {
1295  TVector3 vect1(part1->Px(), part1->Py(), part1->Pz());
1296  TVector3 vect2(part2->Px(), part2->Py(), part2->Pz());
1297  Double_t z = (vect1 * vect2) / (vect2 * vect2);
1298  return z;
1299 }
1300 
1307 Double_t AliAnalysisTaskEmcalLight::GetParallelFraction(const TVector3& vect1, AliVParticle* part2)
1308 {
1309  TVector3 vect2(part2->Px(), part2->Py(), part2->Pz());
1310  Double_t z = (vect1 * vect2) / (vect2 * vect2);
1311  return z;
1312 }
1313 
1322 void AliAnalysisTaskEmcalLight::GetEtaPhiDiff(const AliVTrack *t, const AliVCluster *v, Double_t &phidiff, Double_t &etadiff)
1323 {
1324  phidiff = 999;
1325  etadiff = 999;
1326 
1327  if (!t||!v) return;
1328 
1329  Double_t veta = t->GetTrackEtaOnEMCal();
1330  Double_t vphi = t->GetTrackPhiOnEMCal();
1331 
1332  Float_t pos[3] = {0};
1333  v->GetPosition(pos);
1334  TVector3 cpos(pos);
1335  Double_t ceta = cpos.Eta();
1336  Double_t cphi = cpos.Phi();
1337  etadiff=veta-ceta;
1338  phidiff=TVector2::Phi_mpi_pi(vphi-cphi);
1339 }
1340 
1347 {
1348  Byte_t ret = 0;
1349  if (t->TestBit(BIT(22)) && !t->TestBit(BIT(23)))
1350  ret = 1;
1351  else if (!t->TestBit(BIT(22)) && t->TestBit(BIT(23)))
1352  ret = 2;
1353  else if (t->TestBit(BIT(22)) && t->TestBit(BIT(23)))
1354  ret = 3;
1355  return ret;
1356 }
1357 
1367 Byte_t AliAnalysisTaskEmcalLight::GetTrackType(const AliAODTrack *aodTrack, UInt_t filterBit1, UInt_t filterBit2)
1368 {
1369 
1370  Int_t res = 0;
1371 
1372  if (aodTrack->TestFilterBit(filterBit1)) {
1373  res = 0;
1374  }
1375  else if (aodTrack->TestFilterBit(filterBit2)) {
1376  if ((aodTrack->GetStatus()&AliVTrack::kITSrefit)!=0) {
1377  res = 1;
1378  }
1379  else {
1380  res = 2;
1381  }
1382  }
1383  else {
1384  res = 3;
1385  }
1386 
1387  return res;
1388 }
1389 
1396 {
1397  EBeamType_t b = kpp;
1398  if ((runnumber >= 136833 && runnumber <= 139517) || // LHC10h Run-1 (Pb-Pb)
1399  (runnumber >= 167693 && runnumber <= 170593) || // LHC11h Run-1 (Pb-Pb)
1400  (runnumber >= 244824 && runnumber <= 246994)) { // LHC15o Run-2 (Pb-Pb)
1401  b = kAA;
1402  }
1403  else if ((runnumber > 188356 && runnumber <= 188503) || // LHC12g Run-1 (p-Pb pilot)
1404  (runnumber >= 195164 && runnumber <= 197388) || // LHC13b,c,d,e,f Run-1 (p-Pb)
1405  (runnumber >= 265077 && runnumber <= 267166)) { // LHC16 Run-2 (p-Pb)
1406  b = kpA;
1407  }
1408  return b;
1409 }
1410 
1417 {
1418  if (!fPythiaHeader || !fMCRejectFilter) return kTRUE;
1419 
1420  // Condition 1: Pythia jet / pT-hard > factor
1421  if (fPtHardAndJetPtFactor > 0.) {
1422  AliTLorentzVector jet;
1423 
1424  Int_t nTriggerJets = fPythiaHeader->NTriggerJets();
1425 
1426  AliDebug(1,Form("Njets: %d, pT Hard %f",nTriggerJets, fPtHard));
1427 
1428  Float_t tmpjet[]={0,0,0,0};
1429  for (Int_t ijet = 0; ijet< nTriggerJets; ijet++) {
1430  fPythiaHeader->TriggerJet(ijet, tmpjet);
1431 
1432  jet.SetPxPyPzE(tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3]);
1433 
1434  AliDebug(1,Form("jet %d; pycell jet pT %f",ijet, jet.Pt()));
1435 
1436  //Compare jet pT and pt Hard
1437  if (jet.Pt() > fPtHardAndJetPtFactor * fPtHard) {
1438  AliInfo(Form("Reject jet event with : pT Hard %2.2f, pycell jet pT %2.2f, rejection factor %1.1f\n", fPtHard, jet.Pt(), fPtHardAndJetPtFactor));
1439  return kFALSE;
1440  }
1441  }
1442  }
1443  // end condition 1
1444 
1445  // Condition 2 : Reconstructed EMCal cluster pT / pT-hard > factor
1446  if (fPtHardAndClusterPtFactor > 0.) {
1447  AliClusterContainer* mccluscont = GetClusterContainer(0);
1448  if ((Bool_t)mccluscont) {
1449  for (auto cluster : mccluscont->all()) {// Not cuts applied ; use accept for cuts
1450  Float_t ecluster = cluster->E();
1451 
1452  if (ecluster > (fPtHardAndClusterPtFactor * fPtHard)) {
1453  AliInfo(Form("Reject : ecluster %2.2f, calo %d, factor %2.2f, ptHard %f",ecluster,cluster->GetType(),fPtHardAndClusterPtFactor,fPtHard));
1454  return kFALSE;
1455  }
1456  }
1457  }
1458  }
1459  // end condition 2
1460 
1461  // condition 3 : Reconstructed track pT / pT-hard >factor
1462  if (fPtHardAndTrackPtFactor > 0.) {
1463  AliMCParticleContainer* mcpartcont = dynamic_cast<AliMCParticleContainer*>(GetParticleContainer(0));
1464  if ((Bool_t)mcpartcont) {
1465  for (auto mctrack : mcpartcont->all()) {// Not cuts applied ; use accept for cuts
1466  Float_t trackpt = mctrack->Pt();
1467  if (trackpt > (fPtHardAndTrackPtFactor * fPtHard) ) {
1468  AliInfo(Form("Reject : track %2.2f, factor %2.2f, ptHard %f", trackpt, fPtHardAndTrackPtFactor, fPtHard));
1469  return kFALSE;
1470  }
1471  }
1472  }
1473  }
1474  // end condition 3
1475 
1476  return kTRUE;
1477 }
1478 
1488 {
1489  Double_t dphi = -999;
1490  const Double_t tpi = TMath::TwoPi();
1491 
1492  if (phia < 0) phia += tpi;
1493  else if (phia > tpi) phia -= tpi;
1494  if (phib < 0) phib += tpi;
1495  else if (phib > tpi) phib -= tpi;
1496  dphi = phib - phia;
1497  if (dphi < rangeMin) dphi += tpi;
1498  else if (dphi > rangeMax) dphi -= tpi;
1499 
1500  return dphi;
1501 }
1502 
1512 void AliAnalysisTaskEmcalLight::GenerateFixedBinArray(int n, double min, double max, std::vector<double>& array, bool last)
1513 {
1514  double binWidth = (max - min) / n;
1515  double v = min;
1516  if (last) n++;
1517  for (int i = 0; i < n; i++) {
1518  array.push_back(v);
1519  v += binWidth;
1520  }
1521 }
1522 
1532 std::vector<double> AliAnalysisTaskEmcalLight::GenerateFixedBinArray(int n, double min, double max, bool last)
1533 {
1534  std::vector<double> array;
1535  GenerateFixedBinArray(n, min, max, array, last);
1536  return array;
1537 }
1538 
1548 void AliAnalysisTaskEmcalLight::GenerateLogFixedBinArray(int n, double min, double max, std::vector<double>& array, bool last)
1549 {
1550  if (min <= 0 || max < min) {
1551  AliErrorClassStream() << "Cannot generate a log scale fixed-bin array with limits " << min << ", " << max << std::endl;
1552  return;
1553  }
1554  double binWidth = std::pow(max / min, 1.0 / n);
1555  double v = min;
1556  if (last) n++;
1557  for (int i = 0; i < n; i++) {
1558  array.push_back(v);
1559  v *= binWidth;
1560  }
1561 }
1562 
1572 std::vector<double> AliAnalysisTaskEmcalLight::GenerateLogFixedBinArray(int n, double min, double max, bool last)
1573 {
1574  std::vector<double> array;
1575  GenerateLogFixedBinArray(n, min, max, array, last);
1576  return array;
1577 }
1578 
1579 
1580 TH1* AliAnalysisTaskEmcalLight::GetGeneralTH1(const char* name, bool warn)
1581 {
1582  auto search = fHistograms.find(name);
1583  if (search != fHistograms.end()) {
1584  return search->second;
1585  }
1586  else {
1587  if (warn) AliErrorStream() << "Could not find histogram '" << name << "'" << std::endl;
1588  return nullptr;
1589  }
1590 }
1591 
1592 TH2* AliAnalysisTaskEmcalLight::GetGeneralTH2(const char* name, bool warn)
1593 {
1594  return static_cast<TH2*>(GetGeneralTH1(name, warn));
1595 }
1596 
1597 TProfile* AliAnalysisTaskEmcalLight::GetGeneralTProfile(const char* name, bool warn)
1598 {
1599  return static_cast<TProfile*>(GetGeneralTH1(name, warn));
1600 }
1601 
1603 {
1604  fIsPythia = i;
1605  if (fIsPythia) {
1606  fIsMonteCarlo = kTRUE;
1607  fMCEventHeaderName = "AliGenPythiaEventHeader";
1608  }
1609  else {
1610  if (fMCEventHeaderName == "AliGenPythiaEventHeader") {
1611  fMCEventHeaderName = "";
1612  }
1613  }
1614 }
1615 
1617 {
1618  TClass gen_header_class(name);
1619  if (gen_header_class.InheritsFrom("AliGenEventHeader")) {
1620  fMCEventHeaderName = name;
1621  }
1622  else {
1623  AliWarningStream() << "Class name '" << name << "' does not inherit from 'AliGenEventHeader'. Not setting it." << std::endl;
1624  }
1625 }
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)