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