AliPhysics  9c66e61 (9c66e61)
AliAnalysisTaskEmcalLight.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2016, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 #include <sstream>
16 #include <array>
17 #include <memory>
18 
19 #include <RVersion.h>
20 #include <TClonesArray.h>
21 #include <TList.h>
22 #include <TObject.h>
23 #include <TH1F.h>
24 #include <TH2F.h>
25 #include <TProfile.h>
26 #include <TSystem.h>
27 #include <TFile.h>
28 #include <TChain.h>
29 #include <TKey.h>
30 
31 #include "AliGenCocktailEventHeader.h"
32 #include "AliStack.h"
33 #include "AliAODEvent.h"
34 #include "AliAnalysisManager.h"
35 #include "AliCentrality.h"
36 #include "AliEMCALGeometry.h"
37 #include "AliESDEvent.h"
38 #include "AliEmcalParticle.h"
39 #include "AliEventplane.h"
40 #include "AliInputEventHandler.h"
41 #include "AliLog.h"
42 #include "AliMCParticle.h"
43 #include "AliVCluster.h"
44 #include "AliVEventHandler.h"
45 #include "AliVParticle.h"
46 #include "AliAODTrack.h"
47 #include "AliVCaloTrigger.h"
48 #include "AliGenPythiaEventHeader.h"
49 #include "AliGenEventHeader.h"
50 #include "AliAODMCHeader.h"
51 #include "AliMCEvent.h"
52 #include "AliEMCALTriggerPatchInfo.h"
53 
54 #include "AliMultSelection.h"
55 
57 
59 
63 
69  fForceBeamType(kNA),
70  fGeneralHistograms(kFALSE),
71  fCreateHisto(kTRUE),
72  fNeedEmcalGeom(kTRUE),
73  fCentBins(),
74  fCentralityEstimation(kNewCentrality),
75  fIsPythia(kFALSE),
76  fIsMonteCarlo(kFALSE),
77  fMCEventHeaderName(),
78  fCaloCellsName(),
79  fCaloTriggersName(),
80  fCaloTriggerPatchInfoName(),
81  fCentEst("V0M"),
82  fParticleCollArray(),
83  fClusterCollArray(),
84  fTriggerSelectionBitMap(0),
85  fMinCent(-1),
86  fMaxCent(-1),
87  fMinVz(-999),
88  fMaxVz(999),
89  fMaxVzDiff(-1),
90  fMinNVertCont(0),
91  fMinPtHard(-1),
92  fMaxPtHard(-1),
93  fMaxMinimumBiasPtHard(-1),
94  fAcceptedTriggerClasses(),
95  fRejectedTriggerClasses(),
96  fMCRejectFilter(kFALSE),
97  fPtHardAndJetPtFactor(0.),
98  fPtHardAndClusterPtFactor(0.),
99  fPtHardAndTrackPtFactor(0.),
100  fSwitchOffLHC15oFaultyBranches(kFALSE),
101  fEventSelectionAfterRun(kFALSE),
102  fSelectGeneratorName(),
103  fMinimumEventWeight(1e-6),
104  fMaximumEventWeight(1e6),
105  fInhibit(kFALSE),
106  fLocalInitialized(kFALSE),
107  fDataType(kAOD),
108  fGeom(0),
109  fCaloCells(0),
110  fCaloTriggers(0),
111  fTriggerPatchInfo(0),
112  fCent(-1),
113  fCentBin(-1),
114  fEPV0(-1.0),
115  fEPV0A(-1.0),
116  fEPV0C(-1.0),
117  fNVertCont(0),
118  fNVertSPDCont(0),
119  fFiredTriggerBitMap(0),
120  fFiredTriggerClasses(),
121  fBeamType(kNA),
122  fMCHeader(0),
123  fPythiaHeader(0),
124  fUseXsecFromHeader(false),
125  fPtHardBin(0),
126  fPtHard(0),
127  fNTrials(0),
128  fXsection(0),
129  fEventWeight(1),
130  fGeneratorName(),
131  fOutput(0),
132  fHistograms()
133 {
134  fVertex[0] = 0;
135  fVertex[1] = 0;
136  fVertex[2] = 0;
137  fVertexSPD[0] = 0;
138  fVertexSPD[1] = 0;
139  fVertexSPD[2] = 0;
140 }
141 
153  AliAnalysisTaskSE(name),
154  fForceBeamType(kNA),
155  fGeneralHistograms(kFALSE),
156  fCreateHisto(kTRUE),
157  fNeedEmcalGeom(kTRUE),
158  fCentBins(6),
159  fCentralityEstimation(kNewCentrality),
160  fIsPythia(kFALSE),
161  fIsMonteCarlo(kFALSE),
162  fMCEventHeaderName(),
163  fCaloCellsName(),
164  fCaloTriggersName(),
165  fCaloTriggerPatchInfoName(),
166  fCentEst("V0M"),
167  fParticleCollArray(),
168  fClusterCollArray(),
169  fTriggerSelectionBitMap(0),
170  fMinCent(-1),
171  fMaxCent(-1),
172  fMinVz(-999),
173  fMaxVz(999),
174  fMaxVzDiff(-1),
175  fMinNVertCont(0),
176  fMinPtHard(-1),
177  fMaxPtHard(-1),
178  fMaxMinimumBiasPtHard(-1),
179  fAcceptedTriggerClasses(),
180  fRejectedTriggerClasses(),
181  fMCRejectFilter(kFALSE),
182  fPtHardAndJetPtFactor(0.),
183  fPtHardAndClusterPtFactor(0.),
184  fPtHardAndTrackPtFactor(0.),
185  fSwitchOffLHC15oFaultyBranches(kFALSE),
186  fEventSelectionAfterRun(kFALSE),
187  fSelectGeneratorName(),
188  fMinimumEventWeight(1e-6),
189  fMaximumEventWeight(1e6),
190  fInhibit(kFALSE),
191  fLocalInitialized(kFALSE),
192  fDataType(kAOD),
193  fGeom(0),
194  fCaloCells(0),
195  fCaloTriggers(0),
196  fTriggerPatchInfo(0),
197  fCent(0),
198  fCentBin(-1),
199  fEPV0(-1.0),
200  fEPV0A(-1.0),
201  fEPV0C(-1.0),
202  fNVertCont(0),
203  fNVertSPDCont(0),
204  fFiredTriggerBitMap(0),
205  fFiredTriggerClasses(),
206  fBeamType(kNA),
207  fMCHeader(0),
208  fPythiaHeader(0),
209  fUseXsecFromHeader(false),
210  fPtHardBin(0),
211  fPtHard(0),
212  fNTrials(0),
213  fXsection(0),
214  fEventWeight(1),
215  fGeneratorName(),
216  fOutput(0),
217  fHistograms()
218 {
219  fVertex[0] = 0;
220  fVertex[1] = 0;
221  fVertex[2] = 0;
222  fVertexSPD[0] = 0;
223  fVertexSPD[1] = 0;
224  fVertexSPD[2] = 0;
225 
226  fCentBins[0] = 0;
227  fCentBins[1] = 10;
228  fCentBins[2] = 30;
229  fCentBins[3] = 50;
230  fCentBins[4] = 90;
231  fCentBins[5] = 100;
232 
233  if (fCreateHisto) DefineOutput(1, TList::Class());
234 }
235 
240 {
241  for (auto cont_it : fParticleCollArray) delete cont_it.second;
242  for (auto cont_it : fClusterCollArray) delete cont_it.second;
243 }
244 
265 {
266  if (fInhibit) {
267  AliWarningStream() << "The execution of this task is inhibited. Returning." << std::endl;
268  return;
269  }
270  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
271  if (mgr) {
272  AliVEventHandler *evhand = mgr->GetInputEventHandler();
273  if (evhand) {
274  if (evhand->InheritsFrom("AliESDInputHandler")) {
275  fDataType = kESD;
276  }
277  else {
278  fDataType = kAOD;
279  }
280  }
281  else {
282  AliError("Event handler not found!");
283  }
284  }
285  else {
286  AliError("Analysis manager not found!");
287  }
288 
289  if (!fCreateHisto)
290  return;
291 
292  OpenFile(1);
293  fOutput = new TList();
294  fOutput->SetOwner(); // @suppress("Ambiguous problem")
295 
297 
298  if (!fGeneralHistograms) return;
299 
300  TH1* h = nullptr;
301 
302  if (fIsMonteCarlo) {
303  auto weight_bins = GenerateLogFixedBinArray(1000, fMinimumEventWeight, fMaximumEventWeight, true);
304 
305  h = new TH1F("fHistEventsVsPtHard", "fHistEventsVsPtHard", 1000, 0, 1000);
306  h->GetXaxis()->SetTitle("#it{p}_{T,hard} (GeV/#it{c})");
307  h->GetYaxis()->SetTitle("events");
308  fOutput->Add(h);
309  fHistograms["fHistEventsVsPtHard"] = h;
310 
311  h = new TH1F("fHistTrialsVsPtHard", "fHistTrialsVsPtHard", 1000, 0, 1000);
312  h->GetXaxis()->SetTitle("#it{p}_{T,hard} (GeV/#it{c})");
313  h->GetYaxis()->SetTitle("trials");
314  fOutput->Add(h);
315  fHistograms["fHistTrialsVsPtHard"] = h;
316 
317  h = new TProfile("fHistXsection", "fHistXsection", 50, 0, 50);
318  h->GetXaxis()->SetTitle("#it{p}_{T,hard} bin");
319  h->GetYaxis()->SetTitle("total integrated cross section (mb)");
320  fOutput->Add(h);
321  fHistograms["fHistXsection"] = h;
322 
323  h = new TH1F("fHistXsectionDistribution", "fHistXsectionDistribution", 1000, &weight_bins[0]);
324  h->GetXaxis()->SetTitle("total integrated cross section (mb)");
325  h->GetYaxis()->SetTitle("events");
326  fOutput->Add(h);
327  fHistograms["fHistXsectionDistribution"] = h;
328 
329  h = new TH1F("fHistEventWeights", "fHistEventWeights", 1000, &weight_bins[0]);
330  h->GetXaxis()->SetTitle("weight");
331  h->GetYaxis()->SetTitle("events");
332  fOutput->Add(h);
333  fHistograms["fHistEventWeights"] = h;
334 
335  h = new TH2F("fHistEventWeightsVsPtHard", "fHistEventWeightsVsPtHard", 1000, 0, 1000, 1000, &weight_bins[0]);
336  h->GetXaxis()->SetTitle("#it{p}_{T,hard} (GeV/#it{c})");
337  h->GetYaxis()->SetTitle("event weight");
338  fOutput->Add(h);
339  fHistograms["fHistEventWeightsVsPtHard"] = h;
340 
341  h = new TH1F("fHistEventsVsPtHardNoSel", "fHistEventsVsPtHardNoSel", 1000, 0, 1000);
342  h->GetXaxis()->SetTitle("#it{p}_{T,hard} (GeV/#it{c})");
343  h->GetYaxis()->SetTitle("events");
344  fOutput->Add(h);
345  fHistograms["fHistEventsVsPtHardNoSel"] = h;
346 
347  h = new TH1F("fHistTrialsVsPtHardNoSel", "fHistTrialsVsPtHardNoSel", 1000, 0, 1000);
348  h->GetXaxis()->SetTitle("#it{p}_{T,hard} (GeV/#it{c})");
349  h->GetYaxis()->SetTitle("trials");
350  fOutput->Add(h);
351  fHistograms["fHistTrialsVsPtHardNoSel"] = h;
352 
353  h = new TProfile("fHistXsectionNoSel", "fHistXsectionNoSel", 50, 0, 50);
354  h->GetXaxis()->SetTitle("#it{p}_{T,hard} bin");
355  h->GetYaxis()->SetTitle("total integrated cross section (mb)");
356  fOutput->Add(h);
357  fHistograms["fHistXsectionNoSel"] = h;
358 
359  h = new TH1F("fHistXsectionDistributionNoSel", "fHistXsectionDistributionNoSel", 1000, &weight_bins[0]);
360  h->GetXaxis()->SetTitle("total integrated cross section (mb)");
361  h->GetYaxis()->SetTitle("events");
362  fOutput->Add(h);
363  fHistograms["fHistXsectionDistributionNoSel"] = h;
364 
365  h = new TH1F("fHistEventWeightsNoSel", "fHistEventWeightsNoSel", 1000, &weight_bins[0]);
366  h->GetXaxis()->SetTitle("weight");
367  h->GetYaxis()->SetTitle("events");
368  fOutput->Add(h);
369  fHistograms["fHistEventWeightsNoSel"] = h;
370 
371  h = new TH2F("fHistEventWeightsVsPtHardNoSel", "fHistEventWeightsVsPtHardNoSel", 1000, 0, 1000, 1000, &weight_bins[0]);
372  h->GetXaxis()->SetTitle("#it{p}_{T,hard} (GeV/#it{c})");
373  h->GetYaxis()->SetTitle("event weight");
374  fOutput->Add(h);
375  fHistograms["fHistEventWeightsVsPtHardNoSel"] = h;
376 
377  h = new TH1F("fHistTrialsExternalFile", "fHistTrialsExternalFile", 50, 0, 50);
378  h->GetXaxis()->SetTitle("#it{p}_{T,hard} bin");
379  h->GetYaxis()->SetTitle("trials");
380  fOutput->Add(h);
381  fHistograms["fHistTrialsExternalFile"] = h;
382 
383  h = new TH1F("fHistEventsExternalFile", "fHistEventsExternalFile", 50, 0, 50);
384  h->GetXaxis()->SetTitle("#it{p}_{T,hard} bin");
385  h->GetYaxis()->SetTitle("total events");
386  fOutput->Add(h);
387  fHistograms["fHistEventsExternalFile"] = h;
388 
389  h = new TProfile("fHistXsectionExternalFile", "fHistXsectionExternalFile", 50, 0, 50);
390  h->GetXaxis()->SetTitle("#it{p}_{T,hard} bin");
391  h->GetYaxis()->SetTitle("total integrated cross section (mb)");
392  fOutput->Add(h);
393  fHistograms["fHistXsectionExternalFile"] = h;
394  }
395 
396  h = new TH1F("fHistZVertex","Z vertex position", 60, -30, 30);
397  h->GetXaxis()->SetTitle("V_{#it{z}}");
398  h->GetYaxis()->SetTitle("counts");
399  fOutput->Add(h);
400  fHistograms["fHistZVertex"] = h;
401 
402  h = new TH1F("fHistZVertexNoSel","Z vertex position (no event selection)", 60, -30, 30);
403  h->GetXaxis()->SetTitle("V_{#it{z}}");
404  h->GetYaxis()->SetTitle("counts");
405  fOutput->Add(h);
406  fHistograms["fHistZVertexNoSel"] = h;
407 
409  h = new TH1F("fHistCentrality","Event centrality distribution", 100, 0, 100);
410  h->GetXaxis()->SetTitle("Centrality (%)");
411  h->GetYaxis()->SetTitle("counts");
412  fOutput->Add(h);
413  fHistograms["fHistCentrality"] = h;
414 
415  h = new TH1F("fHistCentralityNoSel","Event centrality distribution (no event selection)", 100, 0, 100);
416  h->GetXaxis()->SetTitle("Centrality (%)");
417  h->GetYaxis()->SetTitle("counts");
418  fOutput->Add(h);
419  fHistograms["fHistCentralityNoSel"] = h;
420  }
421 
422  if (fForceBeamType != kpp) {
423  h = new TH1F("fHistEventPlane","Event plane", 120, -TMath::Pi(), TMath::Pi());
424  h->GetXaxis()->SetTitle("event plane");
425  h->GetYaxis()->SetTitle("counts");
426  fOutput->Add(h);
427  fHistograms["fHistEventPlane"] = h;
428 
429  h = new TH1F("fHistEventPlaneNoSel","Event plane (no event selection)", 120, -TMath::Pi(), TMath::Pi());
430  h->GetXaxis()->SetTitle("event plane");
431  h->GetYaxis()->SetTitle("counts");
432  fOutput->Add(h);
433  fHistograms["fHistEventPlaneNoSel"] = h;
434  }
435 
436  h = new TH1F("fHistEventRejection","Reasons to reject event",30,0,30);
437 #if ROOT_VERSION_CODE < ROOT_VERSION(6,4,2)
438  h->SetBit(TH1::kCanRebin);
439 #else
440  h->SetCanExtend(TH1::kAllAxes);
441 #endif
442  std::array<std::string, 10> labels = {"PhysSel", "Evt Gen Name", "Trg class (acc)", "Trg class (rej)", "Cent", "vertex contr.", "Vz", "VzSPD", "SelPtHardBin", "MCOutlier"};
443  int i = 1;
444  for (auto label : labels) {
445  h->GetXaxis()->SetBinLabel(i, label.c_str());
446  i++;
447  }
448  h->GetYaxis()->SetTitle("counts");
449  fOutput->Add(h);
450  fHistograms["fHistEventRejection"] = h;
451 
452  h = new TH1F("fHistTriggerClasses","fHistTriggerClasses",3,0,3);
453 #if ROOT_VERSION_CODE < ROOT_VERSION(6,4,2)
454  h->SetBit(TH1::kCanRebin);
455 #else
456  h->SetCanExtend(TH1::kAllAxes);
457 #endif
458  fOutput->Add(h);
459  fHistograms["fHistTriggerClasses"] = h;
460 
461  h = new TH1F("fHistTriggerClassesNoSel","fHistTriggerClassesNoSel",3,0,3);
462 #if ROOT_VERSION_CODE < ROOT_VERSION(6,4,2)
463  h->SetBit(TH1::kCanRebin);
464 #else
465  h->SetCanExtend(TH1::kAllAxes);
466 #endif
467  fOutput->Add(h);
468  fHistograms["fHistTriggerClassesNoSel"] = h;
469 
470  h = new TH1F("fHistEventCount","fHistEventCount",2,0,2);
471  h->GetXaxis()->SetBinLabel(1,"Accepted");
472  h->GetXaxis()->SetBinLabel(2,"Rejected");
473  h->GetYaxis()->SetTitle("counts");
474  fOutput->Add(h);
475  fHistograms["fHistEventCount"] = h;
476 
477  PostData(1, fOutput);
478 }
479 
495 {
496  if (eventSelected) {
497  if (fIsMonteCarlo) {
498  GetGeneralTH1("fHistEventsVsPtHard", true)->Fill(fPtHard);
499  GetGeneralTH1("fHistTrialsVsPtHard", true)->Fill(fPtHard, fNTrials);
500  GetGeneralTH1("fHistEventWeights", true)->Fill(fEventWeight);
501  GetGeneralTH2("fHistEventWeightsVsPtHard", true)->Fill(fPtHard, fEventWeight);
502  GetGeneralTH1("fHistXsectionDistribution", true)->Fill(fXsection);
503  GetGeneralTProfile("fHistXsection", true)->Fill(fPtHardBin, fXsection);
504  }
505 
506  GetGeneralTH1("fHistZVertex")->Fill(fVertex[2]);
507 
508  TH1* hCent = GetGeneralTH1("fHistCentrality");
509  if (hCent) hCent->Fill(fCent);
510 
511  TH1* hEventPlane = GetGeneralTH1("fHistEventPlane");
512  if (hEventPlane) hEventPlane->Fill(fEPV0);
513 
514  TH1* hTriggerClasses = GetGeneralTH1("fHistTriggerClasses");
515  for (auto fired_trg : fFiredTriggerClasses) hTriggerClasses->Fill(fired_trg.c_str(), 1);
516  }
517  else {
518  if (fIsMonteCarlo) {
519  GetGeneralTH1("fHistEventsVsPtHardNoSel", true)->Fill(fPtHard);
520  GetGeneralTH1("fHistTrialsVsPtHardNoSel", true)->Fill(fPtHard, fNTrials);
521  GetGeneralTH1("fHistEventWeightsNoSel", true)->Fill(fEventWeight);
522  GetGeneralTH2("fHistEventWeightsVsPtHardNoSel", true)->Fill(fPtHard, fEventWeight);
523  GetGeneralTH1("fHistXsectionDistributionNoSel", true)->Fill(fXsection);
524  GetGeneralTProfile("fHistXsectionNoSel", true)->Fill(fPtHardBin, fXsection);
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, Bool_t &useXsecFromHeader)
611 {
612 
613  TString file(currFile);
614  xsec = 0;
615  trials = 1;
616 
617  // Determine archive type
618  TString archivetype;
619  std::unique_ptr<TObjArray> walk(file.Tokenize("/"));
620  for(auto t : *walk){
621  TString &tok = static_cast<TObjString *>(t)->String();
622  if(tok.Contains(".zip")){
623  archivetype = tok;
624  Int_t pos = archivetype.Index(".zip");
625  archivetype.Replace(pos, archivetype.Length() - pos, "");
626  }
627  }
628  if(archivetype.Length()){
629  AliDebugStream(1) << "Auto-detected archive type " << archivetype << std::endl;
630  Ssiz_t pos1 = file.Index(archivetype,archivetype.Length(),0,TString::kExact);
631  Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
632  Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
633  file.Replace(pos+1,pos2-pos1,"");
634  } else {
635  // not an archive take the basename....
636  file.ReplaceAll(gSystem->BaseName(file.Data()),"");
637  }
638  AliDebugStream(1) << "File name: " << file << std::endl;
639 
640  // Build virtual file name
641  // Support for train tests
642  TString virtualFileName;
643  if(file.Contains("__alice")){
644  TString tmp(file);
645  Int_t pos = tmp.Index("__alice");
646  tmp.Replace(0, pos, "");
647  tmp.ReplaceAll("__", "/");
648  // cut out tag for archive and root file
649  // this needs a determin
650  std::unique_ptr<TObjArray> toks(tmp.Tokenize("/"));
651  TString tag = "_" + archivetype;
652  for(auto t : *toks){
653  TString &path = static_cast<TObjString *>(t)->String();
654  if(path.Contains(tag)){
655  Int_t posTag = path.Index(tag);
656  path.Replace(posTag, path.Length() - posTag, "");
657  }
658  virtualFileName += "/" + path;
659  }
660  } else {
661  virtualFileName = file;
662  }
663 
664  AliDebugStream(1) << "Physical file name " << file << ", virtual file name " << virtualFileName << std::endl;
665 
666  // Get the pt hard bin
667  TString strPthard(virtualFileName);
668 
669  /*
670  // Dead code - to be removed after testing phase
671  // Procedure will fail for everything else than the expected path name
672  strPthard.Remove(strPthard.Last('/'));
673  strPthard.Remove(strPthard.Last('/'));
674  if (strPthard.Contains("AOD")) strPthard.Remove(strPthard.Last('/'));
675  strPthard.Remove(0,strPthard.Last('/')+1);
676  if (strPthard.IsDec()) pthard = strPthard.Atoi();
677  else
678  AliWarningStream() << "Could not extract file number from path " << strPthard << std::endl;
679  */
680 
681  // New implementation : pattern matching
682  // Reason: Implementation valid only for old productions (new productions swap run number and pt-hard bin)
683  // Idea: Don't use the position in the string but the match different informations
684  // + Year clearly 2000+
685  // + Run number can be match to the one in the event
686  // + If we know it is not year or run number, it must be the pt-hard bin if we start from the beginning
687  // The procedure is only valid for the current implementations and unable to detect non-pt-hard bins
688  // It will also fail in case of arbitrary file names
689 
690  bool binfound = false;
691  std::unique_ptr<TObjArray> tokens(strPthard.Tokenize("/"));
692  for(auto t : *tokens) {
693  TString &tok = static_cast<TObjString *>(t)->String();
694  if(tok.IsDec()){
695  Int_t number = tok.Atoi();
696  if(number > 2000 && number < 3000){
697  // Year
698  continue;
699  } else if(number == fInputHandler->GetEvent()->GetRunNumber()){
700  // Run number
701  continue;
702  } else {
703  if(!binfound){
704  // the first number that is not one of the two must be the pt-hard bin
705  binfound = true;
706  pthard = number;
707  break;
708  }
709  }
710  }
711  }
712  if(!binfound) {
713  AliErrorStream() << "Could not extract file number from path " << strPthard << std::endl;
714  } else {
715  AliInfoStream() << "Auto-detecting pt-hard bin " << pthard << std::endl;
716  }
717 
718  AliInfoStream() << "File: " << file << std::endl;
719 
720  // 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
721  std::unique_ptr<TFile> fxsec(TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")));
722 
723  if (!fxsec) {
724  // next trial fetch the histgram file
725  fxsec = std::unique_ptr<TFile>(TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root")));
726  if (!fxsec){
727  AliErrorStream() << "Failed reading cross section from file " << file << std::endl;
728  useXsecFromHeader = true;
729  return kFALSE; // not a severe condition but inciate that we have no information
730  }
731  else {
732  // find the tlist we want to be independtent of the name so use the Tkey
733  TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
734  if (!key) return kFALSE;
735  TList *list = dynamic_cast<TList*>(key->ReadObj());
736  if (!list) return kFALSE;
737  TProfile *xSecHist = static_cast<TProfile*>(list->FindObject("h1Xsec"));
738  // check for failure
739  if(!xSecHist->GetEntries()) {
740  // No cross seciton information available - fall back to raw
741  AliErrorStream() << "No cross section information available in file " << fxsec->GetName() <<" - fall back to cross section in PYTHIA header" << std::endl;
742  useXsecFromHeader = true;
743  } else {
744  // Cross section histogram filled - take it from there
745  xsec = xSecHist->GetBinContent(1);
746  if(!xsec) AliErrorStream() << GetName() << ": Cross section 0 for file " << file << std::endl;
747  useXsecFromHeader = false;
748  }
749  trials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
750  }
751  } else { // no tree pyxsec.root
752  TTree *xtree = (TTree*)fxsec->Get("Xsection");
753  if (!xtree) return kFALSE;
754  UInt_t ntrials = 0;
755  Double_t xsection = 0;
756  xtree->SetBranchAddress("xsection",&xsection);
757  xtree->SetBranchAddress("ntrials",&ntrials);
758  xtree->GetEntry(0);
759  trials = ntrials;
760  xsec = xsection;
761  }
762  return kTRUE;
763 }
764 
779 {
780  if (!fIsPythia) return kTRUE;
781 
782  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
783  if (!tree) {
784  AliError(Form("%s - UserNotify: No current tree!",GetName()));
785  return kFALSE;
786  }
787 
788  Float_t xsection = 0;
789  Float_t trials = 0;
790  Int_t pthardbin = 0;
791  Bool_t useXsecFromHeader = false;
792 
793  TFile *curfile = tree->GetCurrentFile();
794  if (!curfile) {
795  AliError(Form("%s - UserNotify: No current file!",GetName()));
796  return kFALSE;
797  }
798 
799  TChain *chain = dynamic_cast<TChain*>(tree);
800  if (chain) tree = chain->GetTree();
801 
802  Int_t nevents = tree->GetEntriesFast();
803 
804  Bool_t res = PythiaInfoFromFile(curfile->GetName(), xsection, trials, pthardbin, useXsecFromHeader);
805 
806  fPtHardBin = pthardbin >= 0 ? pthardbin : 0;
807  fUseXsecFromHeader = useXsecFromHeader;
808 
809  if (!res) return kTRUE;
810 
812  GetGeneralTH1("fHistTrialsExternalFile", true)->Fill(fPtHardBin, trials);
813  if(!useXsecFromHeader) GetGeneralTProfile("fHistXsectionExternalFile", true)->Fill(fPtHardBin, xsection);
814  GetGeneralTH1("fHistEventsExternalFile", true)->Fill(fPtHardBin, nevents);
815  }
816 
817  return kTRUE;
818 }
819 
831 {
832  if (!InputEvent()) {
833  AliError(Form("%s: Could not retrieve event! Returning!", GetName()));
834  return;
835  }
836 
837  if (fNeedEmcalGeom && !fGeom) {
838  fGeom = AliEMCALGeometry::GetInstanceFromRunNumber(InputEvent()->GetRunNumber());
839  if (!fGeom) {
840  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()));
841  return;
842  }
843  }
844 
845  if (fSwitchOffLHC15oFaultyBranches && dynamic_cast<AliAODEvent*>(InputEvent())) {
846  TTree *aodTree = AliAnalysisManager::GetAnalysisManager()->GetTree();
847  aodTree->SetBranchStatus("D0toKpi.fPx", 0);
848  aodTree->SetBranchStatus("D0toKpi.fPy", 0);
849  aodTree->SetBranchStatus("D0toKpi.fPz", 0);
850  aodTree->SetBranchStatus("D0toKpi.fd0", 0);
851  aodTree->SetBranchStatus("Charm3Prong.fPx", 0);
852  aodTree->SetBranchStatus("Charm3Prong.fPy", 0);
853  aodTree->SetBranchStatus("Charm3Prong.fPz", 0);
854  aodTree->SetBranchStatus("Charm3Prong.fd0", 0);
855  aodTree->SetBranchStatus("Dstar.fPx", 0);
856  aodTree->SetBranchStatus("Dstar.fPy", 0);
857  aodTree->SetBranchStatus("Dstar.fPz", 0);
858  aodTree->SetBranchStatus("Dstar.fd0", 0);
859  }
860 
861  //Load all requested track branches - each container knows name already
862  for (auto cont_it : fParticleCollArray) cont_it.second->SetArray(InputEvent());
863 
864  //Load all requested cluster branches - each container knows name already
865  for (auto cont_it : fClusterCollArray) cont_it.second->SetArray(InputEvent());
866 
867  if (!fCaloCellsName.IsNull() && !fCaloCells) {
868  fCaloCells = dynamic_cast<AliVCaloCells*>(InputEvent()->FindListObject(fCaloCellsName));
869  if (!fCaloCells) {
870  AliError(Form("%s: Could not retrieve cells %s!", GetName(), fCaloCellsName.Data()));
871  return;
872  }
873  }
874 
875  if (!fCaloTriggersName.IsNull() && !fCaloTriggers) {
876  fCaloTriggers = dynamic_cast<AliVCaloTrigger*>(InputEvent()->FindListObject(fCaloTriggersName));
877  if (!fCaloTriggers) {
878  AliError(Form("%s: Could not retrieve calo triggers %s!", GetName(), fCaloTriggersName.Data()));
879  return;
880  }
881  }
882 
883  if (!fCaloTriggerPatchInfoName.IsNull() && !fTriggerPatchInfo) {
884  fTriggerPatchInfo = GetArrayFromEvent(fCaloTriggerPatchInfoName.Data(),"AliEMCALTriggerPatchInfo");
885  if (!fTriggerPatchInfo) {
886  AliError(Form("%s: Could not retrieve calo trigger patch info %s!", GetName(), fCaloTriggerPatchInfoName.Data()));
887  return;
888  }
889 
890  }
891 
892  fLocalInitialized = kTRUE;
893 }
894 
901 {
902  if (fForceBeamType != kNA)
903  return fForceBeamType;
904 
905  AliESDEvent *esd = dynamic_cast<AliESDEvent*>(InputEvent());
906  if (esd) {
907  const AliESDRun *run = esd->GetESDRun();
908  TString beamType = run->GetBeamType();
909  if (beamType == "p-p")
910  return kpp;
911  else if (beamType == "A-A")
912  return kAA;
913  else if (beamType == "p-A")
914  return kpA;
915  else
916  return kNA;
917  } else {
918  Int_t runNumber = InputEvent()->GetRunNumber();
919  // All run number ranges taken from the RCT
920  if ((runNumber >= 136833 && runNumber <= 139517) || // LHC10h
921  (runNumber >= 167693 && runNumber <= 170593) || // LHC11h
922  (runNumber >= 244824 && runNumber <= 246994)) { // LHC15o
923  return kAA;
924  } else if ((runNumber >= 188356 && runNumber <= 188366) || // LHC12g
925  (runNumber >= 195164 && runNumber <= 197388) || // LHC13b-f
926  (runNumber >= 265015 && runNumber <= 267166)) { // LHC16q-t
927  return kpA;
928  } else {
929  return kpp;
930  }
931  }
932 }
933 
956 {
957  TH1* hEventRejection = GetGeneralTH1("fHistEventRejection", true);
958 
960  if (fGeneralHistograms) hEventRejection->Fill("PhysSel",1);
961  return kFALSE;
962  }
963 
964  if (!fSelectGeneratorName.IsNull() && !fGeneratorName.IsNull()) {
965  if (!fGeneratorName.Contains(fSelectGeneratorName)) {
966  if (fGeneralHistograms) hEventRejection->Fill("Evt Gen Name",1);
967  return kFALSE;
968  }
969  }
970 
971  Bool_t acceptedTrgClassFound = kFALSE;
972  if (fAcceptedTriggerClasses.size() > 0) {
973  for (auto acc_trg : fAcceptedTriggerClasses) {
974  std::string teststring(acc_trg);
975  bool fullmatch(false);
976  auto posexact = acc_trg.find("EXACT");
977  if(posexact != std::string::npos) {
978  fullmatch = true;
979  teststring.erase(posexact, 5);
980  }
981  for (auto fired_trg : fFiredTriggerClasses) {
982  bool classmatch = fullmatch ? teststring == fired_trg : fired_trg.find(teststring) != std::string::npos;
983  if (classmatch) {
984  acceptedTrgClassFound = kTRUE;
985  break;
986  }
987  }
988  if (acceptedTrgClassFound) break;
989  }
990 
991  if (!acceptedTrgClassFound) {
992  if (fGeneralHistograms) hEventRejection->Fill("Trg class (acc)",1);
993  return kFALSE;
994  }
995  }
996 
997  if (fRejectedTriggerClasses.size() > 0) {
998  for (auto rej_trg : fRejectedTriggerClasses) {
999  std::string teststring(rej_trg);
1000  bool fullmatch(false);
1001  auto posexact = rej_trg.find("EXACT");
1002  if(posexact != std::string::npos) {
1003  fullmatch = true;
1004  teststring.erase(posexact, 5);
1005  }
1006  for (auto fired_trg : fFiredTriggerClasses) {
1007  bool classmatch = fullmatch ? teststring == fired_trg : fired_trg.find(teststring) != std::string::npos;
1008  if (classmatch) {
1009  if (fGeneralHistograms) hEventRejection->Fill("Trg class (rej)",1);
1010  return kFALSE;
1011  }
1012  }
1013  }
1014  }
1015 
1016  if (fMinCent < fMaxCent && fMaxCent > 0) {
1017  if (fCent < fMinCent || fCent > fMaxCent) {
1018  if (fGeneralHistograms) hEventRejection->Fill("Cent",1);
1019  return kFALSE;
1020  }
1021  }
1022 
1023  if (fNVertCont < fMinNVertCont) {
1024  if (fGeneralHistograms) hEventRejection->Fill("vertex contr.",1);
1025  return kFALSE;
1026  }
1027 
1028  if (fMinVz < fMaxVz) {
1029  if (fVertex[2] < fMinVz || fVertex[2] > fMaxVz) {
1030  if (fGeneralHistograms) hEventRejection->Fill("Vz",1);
1031  return kFALSE;
1032  }
1033  }
1034 
1035  if (fMaxVzDiff >= 0) {
1036  if (fNVertSPDCont > 0) {
1037  Double_t vzSPD = fVertexSPD[2];
1038  Double_t dvertex = TMath::Abs(fVertex[2] - vzSPD);
1039  //if difference larger than fZvertexDiff
1040  if (dvertex > fMaxVzDiff) {
1041  if (fGeneralHistograms) hEventRejection->Fill("VzSPD",1);
1042  return kFALSE;
1043  }
1044  }
1045  }
1046 
1047  if (fMinPtHard >= 0 && fPtHard < fMinPtHard) {
1048  if (fGeneralHistograms) hEventRejection->Fill("SelPtHardBin",1);
1049  return kFALSE;
1050  }
1051 
1052  if (fMaxPtHard >= 0 && fPtHard >= fMaxPtHard) {
1053  if (fGeneralHistograms) hEventRejection->Fill("SelPtHardBin",1);
1054  return kFALSE;
1055  }
1056 
1058  if (fGeneralHistograms) hEventRejection->Fill("SelPtHardBin",1);
1059  return kFALSE;
1060  }
1061 
1062  // Reject filter for MC data
1063  if (!CheckMCOutliers()) {
1064  if (fGeneralHistograms) hEventRejection->Fill("MCOutlier",1);
1065  return kFALSE;
1066  }
1067 
1068  return kTRUE;
1069 }
1070 
1079 TClonesArray *AliAnalysisTaskEmcalLight::GetArrayFromEvent(const char *name, const char *clname)
1080 {
1081  TClonesArray *arr = 0;
1082  TString sname(name);
1083  if (!sname.IsNull()) {
1084  arr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(sname));
1085  if (!arr) {
1086  AliWarning(Form("%s: Could not retrieve array with name %s!", GetName(), name));
1087  return 0;
1088  }
1089  } else {
1090  return 0;
1091  }
1092 
1093  if (!clname)
1094  return arr;
1095 
1096  TString objname(arr->GetClass()->GetName());
1097  TClass cls(objname);
1098  if (!cls.InheritsFrom(clname)) {
1099  AliWarning(Form("%s: Objects of type %s in %s are not inherited from %s!",
1100  GetName(), cls.GetName(), name, clname));
1101  return 0;
1102  }
1103  return arr;
1104 }
1105 
1111 {
1112  fVertex[0] = 0;
1113  fVertex[1] = 0;
1114  fVertex[2] = 0;
1115  fNVertCont = 0;
1116 
1117  fVertexSPD[0] = 0;
1118  fVertexSPD[1] = 0;
1119  fVertexSPD[2] = 0;
1120  fNVertSPDCont = 0;
1121 
1122  fFiredTriggerClasses.clear();
1123  std::stringstream firedClasses(InputEvent()->GetFiredTriggerClasses().Data());
1124  while (firedClasses.good()) {
1125  std::string trgClass;
1126  firedClasses >> trgClass;
1127  if (!trgClass.empty()) fFiredTriggerClasses.push_back(trgClass);
1128  }
1129 
1130  if (fDataType == kESD) {
1131  fFiredTriggerBitMap = static_cast<AliInputEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler())->IsEventSelected();
1132  }
1133  else {
1134  fFiredTriggerBitMap = static_cast<AliVAODHeader*>(InputEvent()->GetHeader())->GetOfflineTrigger();
1135  }
1136 
1137  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
1138  if (vert) {
1139  vert->GetXYZ(fVertex);
1140  fNVertCont = vert->GetNContributors();
1141  }
1142 
1143  const AliVVertex *vertSPD = InputEvent()->GetPrimaryVertexSPD();
1144  if (vertSPD) {
1145  vertSPD->GetXYZ(fVertexSPD);
1146  fNVertSPDCont = vertSPD->GetNContributors();
1147  }
1148 
1149  fBeamType = GetBeamType();
1150 
1151  fCent = 99;
1152  fCentBin = -1;
1153  fEPV0 = -999;
1154  fEPV0A = -999;
1155  fEPV0C = -999;
1156 
1158  // New centrality estimation (AliMultSelection)
1159  // See https://twiki.cern.ch/twiki/bin/viewauth/ALICE/AliMultSelectionCalibStatus for calibration status period-by-period)
1160  AliMultSelection *MultSelection = static_cast<AliMultSelection*>(InputEvent()->FindListObject("MultSelection"));
1161  if (MultSelection) {
1162  fCent = MultSelection->GetMultiplicityPercentile(fCentEst.Data());
1163  }
1164  else {
1165  AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
1166  }
1167  }
1168  else if (fCentralityEstimation == kOldCentrality) {
1169  // Old centrality estimation (AliCentrality, works only on Run-1 PbPb and pPb)
1170  AliCentrality *aliCent = InputEvent()->GetCentrality();
1171  if (aliCent) {
1172  fCent = aliCent->GetCentralityPercentile(fCentEst.Data());
1173  }
1174  else {
1175  AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
1176  }
1177  }
1178  if (!fCentBins.empty() && fCentralityEstimation != kNoCentrality) {
1179  for (auto cent_it = fCentBins.begin(); cent_it != fCentBins.end() - 1; cent_it++) {
1180  if (fCent >= *cent_it && fCent < *(cent_it+1)) fCentBin = cent_it - fCentBins.begin();
1181  }
1182  }
1183  else {
1184  fCentBin = 0;
1185  }
1186 
1187  if (fBeamType == kAA || fBeamType == kpA ) {
1188  AliEventplane *aliEP = InputEvent()->GetEventplane();
1189  if (aliEP) {
1190  fEPV0 = aliEP->GetEventplane("V0" ,InputEvent());
1191  fEPV0A = aliEP->GetEventplane("V0A",InputEvent());
1192  fEPV0C = aliEP->GetEventplane("V0C",InputEvent());
1193  } else {
1194  AliWarning(Form("%s: Could not retrieve event plane information!", GetName()));
1195  }
1196  }
1197 
1198  if (fIsMonteCarlo && MCEvent()) {
1199  AliGenEventHeader* header = MCEvent()->GenEventHeader();
1200  if (fMCEventHeaderName.IsNull()) {
1201  fMCHeader = header;
1202  }
1203  else {
1204  if (header->InheritsFrom(fMCEventHeaderName)) {
1205  fMCHeader = header;
1206  }
1207  else if (header->InheritsFrom("AliGenCocktailEventHeader")) {
1208  AliGenCocktailEventHeader* cocktailHeader = static_cast<AliGenCocktailEventHeader*>(header);
1209  TList* headers = cocktailHeader->GetHeaders();
1210  for (auto obj : *headers) { // @suppress("Symbol is not resolved")
1211  if (obj->InheritsFrom(fMCEventHeaderName)){
1212  fMCHeader = static_cast<AliGenEventHeader*>(obj);
1213  break;
1214  }
1215  }
1216  }
1217  }
1218  if (fMCHeader) {
1219  fEventWeight = fMCHeader->EventWeight();
1220  if (fIsPythia) {
1221  fPythiaHeader = static_cast<AliGenPythiaEventHeader*>(fMCHeader);
1222  fPtHard = fPythiaHeader->GetPtHard();
1223  fXsection = fPythiaHeader->GetXsection();
1224  fNTrials = fPythiaHeader->Trials();
1225  if(fUseXsecFromHeader) GetGeneralTProfile("fHistXsectionExternalFile", true)->Fill(fPtHardBin, fXsection);
1226  }
1227  }
1228  }
1229 
1230  for (auto cont_it : fParticleCollArray) cont_it.second->NextEvent(InputEvent());
1231  for (auto cont_it : fClusterCollArray) cont_it.second->NextEvent(InputEvent());
1232 
1233  return kTRUE;
1234 }
1235 
1244 AliParticleContainer* AliAnalysisTaskEmcalLight::AddParticleContainer(std::string branchName, std::string contName)
1245 {
1246  if (branchName.size() == 0) return 0;
1247 
1248  AliParticleContainer* cont = 0;
1249 
1250  if (branchName == "tracks" || branchName == "Tracks") cont = new AliTrackContainer(branchName.c_str());
1251  else if (branchName == "mcparticles") cont = new AliMCParticleContainer(branchName.c_str());
1252  else cont = new AliParticleContainer(branchName.c_str());
1253 
1254  if (contName.size() > 0) cont->SetName(contName.c_str());
1255 
1256  AdoptParticleContainer(cont);
1257 
1258  return cont;
1259 }
1260 
1269 AliClusterContainer* AliAnalysisTaskEmcalLight::AddClusterContainer(std::string branchName, std::string contName)
1270 {
1271  if (branchName.size() == 0) return 0;
1272 
1273  AliClusterContainer* cont = new AliClusterContainer(branchName.c_str());
1274 
1275  if (contName.size() > 0) cont->SetName(contName.c_str());
1276 
1277  AdoptClusterContainer(cont);
1278 
1279  return cont;
1280 }
1281 
1288 {
1289  std::map<std::string, AliParticleContainer*>::const_iterator cont_it = fParticleCollArray.find(name);
1290  if (cont_it != fParticleCollArray.end()) return cont_it->second;
1291  else return nullptr;
1292 }
1293 
1300 {
1301  std::map<std::string, AliClusterContainer*>::const_iterator cont_it = fClusterCollArray.find(name);
1302  if (cont_it != fClusterCollArray.end()) return cont_it->second;
1303  else return nullptr;
1304 }
1305 
1312 {
1313  if (!(InputEvent()->FindListObject(obj->GetName()))) {
1314  InputEvent()->AddObject(obj);
1315  }
1316  else {
1317  if (!attempt) {
1318  AliFatal(Form("%s: Container with name %s already present. Aborting", GetName(), obj->GetName()));
1319  }
1320  }
1321 }
1322 
1331 {
1332 
1333  if (!fGeom) {
1334  AliWarning(Form("%s - AliAnalysisTaskEmcalBase::IsTrackInEmcalAcceptance - Geometry is not available!", GetName()));
1335  return kFALSE;
1336  }
1337 
1338  Double_t minPhi = fGeom->GetArm1PhiMin() - edges;
1339  Double_t maxPhi = fGeom->GetArm1PhiMax() + edges;
1340 
1341  if (part->Phi() > minPhi && part->Phi() < maxPhi) {
1342  return kTRUE;
1343  }
1344  else {
1345  return kFALSE;
1346  }
1347 }
1348 
1350 {
1351  axis->SetBinLabel(1, "NullObject");
1352  axis->SetBinLabel(2, "Pt");
1353  axis->SetBinLabel(3, "Acceptance");
1354  axis->SetBinLabel(4, "MCLabel");
1355  axis->SetBinLabel(5, "BitMap");
1356  axis->SetBinLabel(6, "HF cut");
1357  axis->SetBinLabel(7, "Bit6");
1358  axis->SetBinLabel(8, "NotHybridTrack");
1359  axis->SetBinLabel(9, "MCFlag");
1360  axis->SetBinLabel(10, "MCGenerator");
1361  axis->SetBinLabel(11, "ChargeCut");
1362  axis->SetBinLabel(12, "MinDistanceTPCSectorEdge");
1363  axis->SetBinLabel(13, "Bit12");
1364  axis->SetBinLabel(14, "IsEMCal");
1365  axis->SetBinLabel(15, "Time");
1366  axis->SetBinLabel(16, "Energy");
1367  axis->SetBinLabel(17, "ExoticCut");
1368  axis->SetBinLabel(18, "Bit17");
1369  axis->SetBinLabel(19, "Area");
1370  axis->SetBinLabel(20, "AreaEmc");
1371  axis->SetBinLabel(21, "ZLeadingCh");
1372  axis->SetBinLabel(22, "ZLeadingEmc");
1373  axis->SetBinLabel(23, "NEF");
1374  axis->SetBinLabel(24, "MinLeadPt");
1375  axis->SetBinLabel(25, "MaxTrackPt");
1376  axis->SetBinLabel(26, "MaxClusterPt");
1377  axis->SetBinLabel(27, "Flavour");
1378  axis->SetBinLabel(28, "TagStatus");
1379  axis->SetBinLabel(29, "MinNConstituents");
1380  axis->SetBinLabel(30, "Bit29");
1381  axis->SetBinLabel(31, "Bit30");
1382  axis->SetBinLabel(32, "Bit31");
1383 }
1384 
1391 Double_t AliAnalysisTaskEmcalLight::GetParallelFraction(AliVParticle* part1, AliVParticle* part2)
1392 {
1393  TVector3 vect1(part1->Px(), part1->Py(), part1->Pz());
1394  TVector3 vect2(part2->Px(), part2->Py(), part2->Pz());
1395  Double_t z = (vect1 * vect2) / (vect2 * vect2);
1396  return z;
1397 }
1398 
1405 Double_t AliAnalysisTaskEmcalLight::GetParallelFraction(const TVector3& vect1, AliVParticle* part2)
1406 {
1407  TVector3 vect2(part2->Px(), part2->Py(), part2->Pz());
1408  Double_t z = (vect1 * vect2) / (vect2 * vect2);
1409  return z;
1410 }
1411 
1420 void AliAnalysisTaskEmcalLight::GetEtaPhiDiff(const AliVTrack *t, const AliVCluster *v, Double_t &phidiff, Double_t &etadiff)
1421 {
1422  phidiff = 999;
1423  etadiff = 999;
1424 
1425  if (!t||!v) return;
1426 
1427  Double_t veta = t->GetTrackEtaOnEMCal();
1428  Double_t vphi = t->GetTrackPhiOnEMCal();
1429 
1430  Float_t pos[3] = {0};
1431  v->GetPosition(pos);
1432  TVector3 cpos(pos);
1433  Double_t ceta = cpos.Eta();
1434  Double_t cphi = cpos.Phi();
1435  etadiff=veta-ceta;
1436  phidiff=TVector2::Phi_mpi_pi(vphi-cphi);
1437 }
1438 
1445 {
1446  Byte_t ret = 0;
1447  if (t->TestBit(BIT(22)) && !t->TestBit(BIT(23)))
1448  ret = 1;
1449  else if (!t->TestBit(BIT(22)) && t->TestBit(BIT(23)))
1450  ret = 2;
1451  else if (t->TestBit(BIT(22)) && t->TestBit(BIT(23)))
1452  ret = 3;
1453  return ret;
1454 }
1455 
1465 Byte_t AliAnalysisTaskEmcalLight::GetTrackType(const AliAODTrack *aodTrack, UInt_t filterBit1, UInt_t filterBit2)
1466 {
1467 
1468  Int_t res = 0;
1469 
1470  if (aodTrack->TestFilterBit(filterBit1)) {
1471  res = 0;
1472  }
1473  else if (aodTrack->TestFilterBit(filterBit2)) {
1474  if ((aodTrack->GetStatus()&AliVTrack::kITSrefit)!=0) {
1475  res = 1;
1476  }
1477  else {
1478  res = 2;
1479  }
1480  }
1481  else {
1482  res = 3;
1483  }
1484 
1485  return res;
1486 }
1487 
1494 {
1495  EBeamType_t b = kpp;
1496  if ((runnumber >= 136833 && runnumber <= 139517) || // LHC10h Run-1 (Pb-Pb)
1497  (runnumber >= 167693 && runnumber <= 170593) || // LHC11h Run-1 (Pb-Pb)
1498  (runnumber >= 244824 && runnumber <= 246994)) { // LHC15o Run-2 (Pb-Pb)
1499  b = kAA;
1500  }
1501  else if ((runnumber > 188356 && runnumber <= 188503) || // LHC12g Run-1 (p-Pb pilot)
1502  (runnumber >= 195164 && runnumber <= 197388) || // LHC13b,c,d,e,f Run-1 (p-Pb)
1503  (runnumber >= 265077 && runnumber <= 267166)) { // LHC16 Run-2 (p-Pb)
1504  b = kpA;
1505  }
1506  return b;
1507 }
1508 
1515 {
1516  if (!fPythiaHeader || !fMCRejectFilter) return kTRUE;
1517 
1518  // Condition 1: Pythia jet / pT-hard > factor
1519  if (fPtHardAndJetPtFactor > 0.) {
1520  AliTLorentzVector jet;
1521 
1522  Int_t nTriggerJets = fPythiaHeader->NTriggerJets();
1523 
1524  AliDebug(1,Form("Njets: %d, pT Hard %f",nTriggerJets, fPtHard));
1525 
1526  Float_t tmpjet[]={0,0,0,0};
1527  for (Int_t ijet = 0; ijet< nTriggerJets; ijet++) {
1528  fPythiaHeader->TriggerJet(ijet, tmpjet);
1529 
1530  jet.SetPxPyPzE(tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3]);
1531 
1532  AliDebug(1,Form("jet %d; pycell jet pT %f",ijet, jet.Pt()));
1533 
1534  //Compare jet pT and pt Hard
1535  if (jet.Pt() > fPtHardAndJetPtFactor * fPtHard) {
1536  AliInfo(Form("Reject jet event with : pT Hard %2.2f, pycell jet pT %2.2f, rejection factor %1.1f\n", fPtHard, jet.Pt(), fPtHardAndJetPtFactor));
1537  return kFALSE;
1538  }
1539  }
1540  }
1541  // end condition 1
1542 
1543  // Condition 2 : Reconstructed EMCal cluster pT / pT-hard > factor
1544  if (fPtHardAndClusterPtFactor > 0.) {
1545  AliClusterContainer* mccluscont = GetClusterContainer(0);
1546  if ((Bool_t)mccluscont) {
1547  for (auto cluster : mccluscont->all()) {// Not cuts applied ; use accept for cuts
1548  Float_t ecluster = cluster->E();
1549 
1550  if (ecluster > (fPtHardAndClusterPtFactor * fPtHard)) {
1551  AliInfo(Form("Reject : ecluster %2.2f, calo %d, factor %2.2f, ptHard %f",ecluster,cluster->GetType(),fPtHardAndClusterPtFactor,fPtHard));
1552  return kFALSE;
1553  }
1554  }
1555  }
1556  }
1557  // end condition 2
1558 
1559  // condition 3 : Reconstructed track pT / pT-hard >factor
1560  if (fPtHardAndTrackPtFactor > 0.) {
1561  AliMCParticleContainer* mcpartcont = dynamic_cast<AliMCParticleContainer*>(GetParticleContainer(0));
1562  if ((Bool_t)mcpartcont) {
1563  for (auto mctrack : mcpartcont->all()) {// Not cuts applied ; use accept for cuts
1564  Float_t trackpt = mctrack->Pt();
1565  if (trackpt > (fPtHardAndTrackPtFactor * fPtHard) ) {
1566  AliInfo(Form("Reject : track %2.2f, factor %2.2f, ptHard %f", trackpt, fPtHardAndTrackPtFactor, fPtHard));
1567  return kFALSE;
1568  }
1569  }
1570  }
1571  }
1572  // end condition 3
1573 
1574  return kTRUE;
1575 }
1576 
1586 {
1587  Double_t dphi = -999;
1588  const Double_t tpi = TMath::TwoPi();
1589 
1590  if (phia < 0) phia += tpi;
1591  else if (phia > tpi) phia -= tpi;
1592  if (phib < 0) phib += tpi;
1593  else if (phib > tpi) phib -= tpi;
1594  dphi = phib - phia;
1595  if (dphi < rangeMin) dphi += tpi;
1596  else if (dphi > rangeMax) dphi -= tpi;
1597 
1598  return dphi;
1599 }
1600 
1610 void AliAnalysisTaskEmcalLight::GenerateFixedBinArray(int n, double min, double max, std::vector<double>& array, bool last)
1611 {
1612  double binWidth = (max - min) / n;
1613  double v = min;
1614  if (last) n++;
1615  for (int i = 0; i < n; i++) {
1616  array.push_back(v);
1617  v += binWidth;
1618  }
1619 }
1620 
1630 std::vector<double> AliAnalysisTaskEmcalLight::GenerateFixedBinArray(int n, double min, double max, bool last)
1631 {
1632  std::vector<double> array;
1633  GenerateFixedBinArray(n, min, max, array, last);
1634  return array;
1635 }
1636 
1646 void AliAnalysisTaskEmcalLight::GenerateLogFixedBinArray(int n, double min, double max, std::vector<double>& array, bool last)
1647 {
1648  if (min <= 0 || max < min) {
1649  AliErrorClassStream() << "Cannot generate a log scale fixed-bin array with limits " << min << ", " << max << std::endl;
1650  return;
1651  }
1652  double binWidth = std::pow(max / min, 1.0 / n);
1653  double v = min;
1654  if (last) n++;
1655  for (int i = 0; i < n; i++) {
1656  array.push_back(v);
1657  v *= binWidth;
1658  }
1659 }
1660 
1670 std::vector<double> AliAnalysisTaskEmcalLight::GenerateLogFixedBinArray(int n, double min, double max, bool last)
1671 {
1672  std::vector<double> array;
1673  GenerateLogFixedBinArray(n, min, max, array, last);
1674  return array;
1675 }
1676 
1677 
1678 TH1* AliAnalysisTaskEmcalLight::GetGeneralTH1(const char* name, bool warn)
1679 {
1680  auto search = fHistograms.find(name);
1681  if (search != fHistograms.end()) {
1682  return search->second;
1683  }
1684  else {
1685  if (warn) AliErrorStream() << "Could not find histogram '" << name << "'" << std::endl;
1686  return nullptr;
1687  }
1688 }
1689 
1690 TH2* AliAnalysisTaskEmcalLight::GetGeneralTH2(const char* name, bool warn)
1691 {
1692  return static_cast<TH2*>(GetGeneralTH1(name, warn));
1693 }
1694 
1695 TProfile* AliAnalysisTaskEmcalLight::GetGeneralTProfile(const char* name, bool warn)
1696 {
1697  return static_cast<TProfile*>(GetGeneralTH1(name, warn));
1698 }
1699 
1701 {
1702  fIsPythia = i;
1703  if (fIsPythia) {
1704  fIsMonteCarlo = kTRUE;
1705  fMCEventHeaderName = "AliGenPythiaEventHeader";
1706  }
1707  else {
1708  if (fMCEventHeaderName == "AliGenPythiaEventHeader") {
1709  fMCEventHeaderName = "";
1710  }
1711  }
1712 }
1713 
1715 {
1716  TClass gen_header_class(name);
1717  if (gen_header_class.InheritsFrom("AliGenEventHeader")) {
1718  fMCEventHeaderName = name;
1719  }
1720  else {
1721  AliWarningStream() << "Class name '" << name << "' does not inherit from 'AliGenEventHeader'. Not setting it." << std::endl;
1722  }
1723 }
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
Bool_t PythiaInfoFromFile(const char *currFile, Float_t &fXsec, Float_t &fTrials, Int_t &pthard, Bool_t &useXsecFromHeader)
int Int_t
Definition: External.C:63
TH2 * GetGeneralTH2(const char *name, bool warn=false)
Double_t fMinPtHard
select minimum pt hard (MC)
unsigned int UInt_t
Definition: External.C:33
TString fMCEventHeaderName
Looks for MC event properties in a particular MC event type (useful for a MC cocktail production) ...
float Float_t
Definition: External.C:68
Base task in the EMCAL framework (lighter version of AliAnalysisTaskEmcal)
Float_t fPtHardAndClusterPtFactor
Factor between ptHard and cluster pT to reject/accept event.
Double_t fMaxVzDiff
upper limit for distance between primary and SPD vertex
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 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_t fUseXsecFromHeader
!Switch for using cross section from header (if not found in pythia file)
bool Bool_t
Definition: External.C:53
std::map< std::string, AliClusterContainer * > fClusterCollArray
cluster collection array
EBeamType_t
Switch for the beam type.
AliVCaloTrigger * fCaloTriggers
!calo triggers
Container structure for EMCAL clusters.
Container for MC-true particles within the EMCAL framework.
TString fGeneratorName
!name of the MC generator used to produce the current event (only AOD)
void AdoptParticleContainer(AliParticleContainer *cont)
virtual Bool_t FillGeneralHistograms(Bool_t eventSelected)
Double_t fMaxVz
max vertex for event selection
Double_t fCent
!event centrality
Double_t fEPV0C
!event plane V0C
Definition: External.C:196
TClonesArray * fTriggerPatchInfo
!trigger patch info array
Double_t fMaxPtHard
select maximum pt hard (MC)
Old centrality estimation (AliCentrality, works only on Run-1 PbPb and pPb)
TList * OpenFile(const char *fname)
Definition: DrawAnaELoss.C:65
Bool_t fNeedEmcalGeom
whether or not the task needs the emcal geometry
static Double_t GetParallelFraction(AliVParticle *part1, AliVParticle *part2)
EDataType_t fDataType
!data type (ESD or AOD)