AliPhysics  c7b8e89 (c7b8e89)
AliAnalysisTaskEmcal.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-1999, 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 <RVersion.h>
16 #include <iostream>
17 #include <memory>
18 
19 #include <TClonesArray.h>
20 #include <AliEmcalList.h>
21 #include <TObject.h>
22 #include <TObjString.h>
23 #include <TH1F.h>
24 #include <TProfile.h>
25 #include <TSystem.h>
26 #include <TFile.h>
27 #include <TChain.h>
28 #include <TKey.h>
29 
30 #include "AliAnalysisTaskEmcal.h"
31 #include "AliAnalysisUtils.h"
32 #include "AliAODEvent.h"
33 #include "AliAODMCHeader.h"
34 #include "AliAODTrack.h"
35 #include "AliAnalysisManager.h"
36 #include "AliCentrality.h"
38 #include "AliEMCALGeometry.h"
39 #include "AliEmcalPythiaInfo.h"
40 #include "AliEMCALTriggerPatchInfo.h"
41 #include "AliESDEvent.h"
42 #include "AliAODInputHandler.h"
43 #include "AliESDInputHandler.h"
44 #include "AliEventCuts.h"
45 #include "AliEventplane.h"
46 #include "AliGenPythiaEventHeader.h"
47 #include "AliGenHerwigEventHeader.h"
48 #include "AliInputEventHandler.h"
49 #include "AliLog.h"
50 #include "AliMCEvent.h"
51 #include "AliMCParticle.h"
52 #include "AliMultiInputEventHandler.h"
53 #include "AliMultSelection.h"
54 #include "AliStack.h"
55 #include "AliVCaloTrigger.h"
56 #include "AliVCluster.h"
57 #include "AliVEventHandler.h"
58 #include "AliVParticle.h"
59 #include "AliNanoAODHeader.h"
61 
63 
65 ClassImp(AliAnalysisTaskEmcal);
67 
69  AliAnalysisTaskSE("AliAnalysisTaskEmcal"),
70  fPythiaInfoName(""),
71  fForceBeamType(kNA),
72  fGeneralHistograms(kFALSE),
73  fLocalInitialized(kFALSE),
74  fFileChanged(kTRUE),
75  fCreateHisto(kTRUE),
76  fCaloCellsName(),
77  fCaloTriggersName(),
78  fCaloTriggerPatchInfoName(),
79  fMinCent(-999),
80  fMaxCent(-999),
81  fMinVz(-999),
82  fMaxVz(999),
83  fMinVertexContrib(1),
84  fTrackPtCut(0),
85  fMinNTrack(0),
86  fZvertexDiff(0.5),
87  fUseAliAnaUtils(kFALSE),
88  fRejectPileup(kFALSE),
89  fTklVsClusSPDCut(kFALSE),
90  fOffTrigger(AliVEvent::kAny),
91  fTrigClass(),
92  fMinBiasRefTrigger("CINT7-B-NOPF-ALLNOTRD"),
93  fTriggerTypeSel(kND),
94  fNbins(250),
95  fMinBinPt(0),
96  fMaxBinPt(250),
97  fMinPtTrackInEmcal(0),
98  fEventPlaneVsEmcal(-1),
99  fMinEventPlane(-1e6),
100  fMaxEventPlane(1e6),
101  fCentEst("V0M"),
102  fRecycleUnusedEmbeddedEventsMode(false),
103  fIsEmbedded(kFALSE),
104  fIsPythia(kFALSE),
105  fIsHerwig(kFALSE),
106  fGetPtHardBinFromName(kTRUE),
107  fSelectPtHardBin(-999),
108  fMinMCLabel(0),
109  fMCLabelShift(0),
110  fNcentBins(4),
111  fNeedEmcalGeom(kTRUE),
112  fParticleCollArray(),
113  fClusterCollArray(),
114  fTriggers(0),
115  fEMCalTriggerMode(kOverlapWithLowThreshold),
116  fUseNewCentralityEstimation(kFALSE),
117  fGeneratePythiaInfoObject(kFALSE),
118  fUsePtHardBinScaling(kFALSE),
119  fUseXsecFromHeader(kFALSE),
120  fMCRejectFilter(kFALSE),
121  fCountDownscaleCorrectedEvents(kFALSE),
122  fUseInternalEventSelection(kFALSE),
123  fPtHardAndJetPtFactor(0.),
124  fPtHardAndClusterPtFactor(0.),
125  fPtHardAndTrackPtFactor(0.),
126  fRunNumber(-1),
127  fAliEventCuts(nullptr),
128  fAliAnalysisUtils(nullptr),
129  fIsEsd(kFALSE),
130  fGeom(nullptr),
131  fTracks(nullptr),
132  fCaloClusters(nullptr),
133  fCaloCells(nullptr),
134  fCaloTriggers(nullptr),
135  fTriggerPatchInfo(nullptr),
136  fCent(0),
137  fCentBin(-1),
138  fEPV0(-1.0),
139  fEPV0A(-1.0),
140  fEPV0C(-1.0),
141  fNVertCont(0),
142  fNVertSPDCont(0),
143  fBeamType(kNA),
144  fPythiaHeader(nullptr),
145  fHerwigHeader(nullptr),
146  fPtHard(0),
147  fPtHardBin(0),
148  fPtHardBinGlobal(-1),
149  fPtHardInitialized(false),
150  fNPtHardBins(11),
151  fPtHardBinning(),
152  fNTrials(0),
153  fXsection(0),
154  fPythiaInfo(nullptr),
155  fOutput(nullptr),
156  fHistEventCount(nullptr),
157  fHistTrialsAfterSel(nullptr),
158  fHistEventsAfterSel(nullptr),
159  fHistXsectionAfterSel(nullptr),
160  fHistTrials(nullptr),
161  fHistEvents(nullptr),
162  fHistXsection(nullptr),
163  fHistPtHard(nullptr),
164  fHistPtHardCorr(nullptr),
165  fHistPtHardCorrGlobal(nullptr),
166  fHistPtHardBinCorr(nullptr),
167  fHistCentrality(nullptr),
168  fHistZVertex(nullptr),
169  fHistEventPlane(nullptr),
170  fHistEventRejection(nullptr),
171  fHistTriggerClasses(nullptr),
172  fHistTriggerClassesCorr(nullptr)
173 {
174  fVertex[0] = 0;
175  fVertex[1] = 0;
176  fVertex[2] = 0;
177  fVertexSPD[0] = 0;
178  fVertexSPD[1] = 0;
179  fVertexSPD[2] = 0;
180 
181  fParticleCollArray.SetOwner(kTRUE);
182  fClusterCollArray.SetOwner(kTRUE);
183 }
184 
186  AliAnalysisTaskSE(name),
187  fPythiaInfoName(""),
189  fGeneralHistograms(kFALSE),
190  fLocalInitialized(kFALSE),
191  fFileChanged(kFALSE),
192  fCreateHisto(histo),
193  fCaloCellsName(),
196  fMinCent(-999),
197  fMaxCent(-999),
198  fMinVz(-999),
199  fMaxVz(999),
201  fTrackPtCut(0),
202  fMinNTrack(0),
203  fZvertexDiff(0.5),
204  fUseAliAnaUtils(kFALSE),
205  fRejectPileup(kFALSE),
206  fTklVsClusSPDCut(kFALSE),
207  fOffTrigger(AliVEvent::kAny),
208  fTrigClass(),
209  fMinBiasRefTrigger("CINT7-B-NOPF-ALLNOTRD"),
211  fNbins(250),
212  fMinBinPt(0),
213  fMaxBinPt(250),
215  fEventPlaneVsEmcal(-1),
216  fMinEventPlane(-1e6),
217  fMaxEventPlane(1e6),
218  fCentEst("V0M"),
220  fIsEmbedded(kFALSE),
221  fIsPythia(kFALSE),
222  fIsHerwig(kFALSE),
223  fGetPtHardBinFromName(kTRUE),
224  fSelectPtHardBin(-999),
225  fMinMCLabel(0),
226  fMCLabelShift(0),
227  fNcentBins(4),
228  fNeedEmcalGeom(kTRUE),
231  fTriggers(0),
235  fUsePtHardBinScaling(kFALSE),
236  fUseXsecFromHeader(kFALSE),
237  fMCRejectFilter(kFALSE),
243  fRunNumber(-1),
246  fIsEsd(kFALSE),
247  fGeom(nullptr),
248  fTracks(nullptr),
253  fCent(0),
254  fCentBin(-1),
255  fEPV0(-1.0),
256  fEPV0A(-1.0),
257  fEPV0C(-1.0),
258  fNVertCont(0),
259  fNVertSPDCont(0),
260  fBeamType(kNA),
263  fPtHard(0),
264  fPtHardBin(0),
265  fPtHardBinGlobal(-1),
266  fPtHardInitialized(false),
267  fNPtHardBins(11),
268  fPtHardBinning(),
269  fNTrials(0),
270  fXsection(0),
271  fPythiaInfo(0),
272  fOutput(nullptr),
290 {
291  fVertex[0] = 0;
292  fVertex[1] = 0;
293  fVertex[2] = 0;
294  fVertexSPD[0] = 0;
295  fVertexSPD[1] = 0;
296  fVertexSPD[2] = 0;
297  fParticleCollArray.SetOwner(kTRUE);
298  fClusterCollArray.SetOwner(kTRUE);
299 
300  if (fCreateHisto) {
301  DefineOutput(1, AliEmcalList::Class());
302  }
303 }
304 
306 {
307  if(fOutput) delete fOutput;
309  if(fPythiaInfo) delete fPythiaInfo;
310 }
311 
313 {
315  if (cont) cont->SetClusPtCut(cut);
316  else AliError(Form("%s in SetClusPtCut(...): container %d not found",GetName(),c));
317 }
318 
320 {
322  if (cont) cont->SetClusTimeCut(min,max);
323  else AliError(Form("%s in SetClusTimeCut(...): container %d not found",GetName(),c));
324 }
325 
327 {
329  if (cont) cont->SetParticlePtCut(cut);
330  else AliError(Form("%s in SetTrackPtCut(...): container %d not found",GetName(),c));
331 
332  fTrackPtCut = cut;
333 }
334 
336 {
338  if (cont) cont->SetParticleEtaLimits(min,max);
339  else AliError(Form("%s in SetTrackPtCut(...): container %d not found",GetName(),c));
340 }
341 
343 {
345  if (cont) cont->SetParticlePhiLimits(min,max);
346  else AliError(Form("%s in SetTrackPhiLimits(...): container %d not found",GetName(),c));
347 }
348 
350 {
351  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
352  if (mgr) {
353  AliVEventHandler *evhand = mgr->GetInputEventHandler();
354  if (evhand) {
355  if (evhand->InheritsFrom("AliESDInputHandler")) {
356  fIsEsd = kTRUE;
357  }
358  else {
359  fIsEsd = kFALSE;
360  }
361  }
362  else {
363  AliError("Event handler not found!");
364  }
365  }
366  else {
367  AliError("Analysis manager not found!");
368  }
369 
370 
371  if (!fCreateHisto)
372  return;
373 
374  OpenFile(1);
375  fOutput = new AliEmcalList();
377  fOutput->SetOwner();
378 
379  if (fForceBeamType == kpp)
380  fNcentBins = 1;
381 
382  if (!fGeneralHistograms)
383  return;
384 
385  if (fIsPythia || fIsHerwig) {
386  fHistTrialsAfterSel = new TH1F("fHistTrialsAfterSel", "fHistTrialsAfterSel", fNPtHardBins, 0, fNPtHardBins);
387  fHistTrialsAfterSel->GetXaxis()->SetTitle("p_{T} hard bin");
388  fHistTrialsAfterSel->GetYaxis()->SetTitle("trials");
390 
391  fHistEventsAfterSel = new TH1F("fHistEventsAfterSel", "fHistEventsAfterSel", fNPtHardBins, 0, fNPtHardBins);
392  fHistEventsAfterSel->GetXaxis()->SetTitle("p_{T} hard bin");
393  fHistEventsAfterSel->GetYaxis()->SetTitle("total events");
395 
396  fHistXsectionAfterSel = new TProfile("fHistXsectionAfterSel", "fHistXsectionAfterSel", fNPtHardBins, 0, fNPtHardBins);
397  fHistXsectionAfterSel->GetXaxis()->SetTitle("p_{T} hard bin");
398  fHistXsectionAfterSel->GetYaxis()->SetTitle("xsection");
400 
401  fHistTrials = new TH1F("fHistTrials", "fHistTrials", fNPtHardBins, 0, fNPtHardBins);
402  fHistTrials->GetXaxis()->SetTitle("p_{T} hard bin");
403  fHistTrials->GetYaxis()->SetTitle("trials");
404  fOutput->Add(fHistTrials);
405 
406  fHistEvents = new TH1F("fHistEvents", "fHistEvents", fNPtHardBins, 0, fNPtHardBins);
407  fHistEvents->GetXaxis()->SetTitle("p_{T} hard bin");
408  fHistEvents->GetYaxis()->SetTitle("total events");
409  fOutput->Add(fHistEvents);
410 
411  fHistXsection = new TProfile("fHistXsection", "fHistXsection", fNPtHardBins, 0, fNPtHardBins);
412  fHistXsection->GetXaxis()->SetTitle("p_{T} hard bin");
413  fHistXsection->GetYaxis()->SetTitle("xsection");
414  fOutput->Add(fHistXsection);
415 
416  // Set the bin labels
417  Bool_t binningAvailable = false;
418  if(fPtHardBinning.GetSize() > 0) {
419  AliInfoStream() << "Using custom pt-hard binning" << std::endl;
420  if(fPtHardBinning.GetSize() == fNPtHardBins + 1) binningAvailable = true;
421  else AliErrorStream() << "Pt-hard binning (" << fPtHardBinning.GetSize() -1 << ") does not match the amount of bins (" << fNPtHardBins << ")" << std::endl;
422  } else {
423  // Check if we fall back to the default binning
424  if(fNPtHardBins == 11) {
425  AliInfoStream() << "11 pt-hard bins - fall back to default binning for bin labels" << std::endl;
426  const Int_t kDefaultPtHardBinning[12] = {0,5,11,21,36,57, 84,117,152,191,234,1000000};
427  fPtHardBinning.Set(12);
428  for(Int_t ib = 0; ib < 12; ib++) fPtHardBinning[ib] = kDefaultPtHardBinning[ib];
429  binningAvailable = true;
430  } else {
431  AliErrorStream() << "No default binning available for " << fNPtHardBins << " pt-hard bins - bin labels will not be set." << std::endl;
432  }
433  }
434 
435  if(binningAvailable){
436  for (Int_t i = 0; i < fNPtHardBins; i++) {
437  fHistTrialsAfterSel->GetXaxis()->SetBinLabel(i+1, Form("%d-%d",fPtHardBinning[i],fPtHardBinning[i+1]));
438  fHistEventsAfterSel->GetXaxis()->SetBinLabel(i+1, Form("%d-%d",fPtHardBinning[i],fPtHardBinning[i+1]));
439  fHistXsectionAfterSel->GetXaxis()->SetBinLabel(i+1, Form("%d-%d",fPtHardBinning[i],fPtHardBinning[i+1]));
440 
441  fHistTrials->GetXaxis()->SetBinLabel(i+1, Form("%d-%d",fPtHardBinning[i],fPtHardBinning[i+1]));
442  fHistXsection->GetXaxis()->SetBinLabel(i+1, Form("%d-%d",fPtHardBinning[i],fPtHardBinning[i+1]));
443  fHistEvents->GetXaxis()->SetBinLabel(i+1, Form("%d-%d",fPtHardBinning[i],fPtHardBinning[i+1]));
444  }
445  } else {
446  AliErrorStream() << "No suitable binning available - skipping bin labels" << std::endl;
447  }
448 
449 
450  fHistPtHard = new TH1F("fHistPtHard", "fHistPtHard", fNbins*2, fMinBinPt, fMaxBinPt*4);
451  fHistPtHard->GetXaxis()->SetTitle("p_{T,hard} (GeV/c)");
452  fHistPtHard->GetYaxis()->SetTitle("counts");
453  fOutput->Add(fHistPtHard);
454 
455  // Pt-hard correlation histograms (debugging)
456  fHistPtHardCorr = new TH2F("fHistPtHardCorr", "Correlation between pt-hard value and binning", fNPtHardBins, -0.5, fNPtHardBins - 0.5, 1000, 0., 500.);
457  fHistPtHardCorr->SetXTitle("p_{T,hard bin}");
458  fHistPtHardCorr->SetYTitle("p_{T,hard value}");
459  fOutput->Add(fHistPtHardCorr);
460 
461  fHistPtHardCorrGlobal = new TH2F("fHistPtHardCorrGlobal", "Correlation between global pt-hard value and binning", fNPtHardBins, -0.5, fNPtHardBins - 0.5, 1000, 0., 500.);
462  fHistPtHardCorrGlobal->SetXTitle("p_{T,hard} bin_{global}");
463  fHistPtHardCorrGlobal->SetYTitle("p_{T,hard} value");
465 
466  fHistPtHardBinCorr = new TH2F("fHistPtHardBinCorr", "Correlation between global and local pt-hard bin", fNPtHardBins, -0.5, fNPtHardBins - 0.5, fNPtHardBins, -0.5, fNPtHardBins);
467  fHistPtHardBinCorr->SetXTitle("p_{T,hard} bin_{local}");
468  fHistPtHardBinCorr->SetYTitle("p_{T,hard} bin_{global}");
470  }
471 
472  fHistZVertex = new TH1F("fHistZVertex","Z vertex position", 60, -30, 30);
473  fHistZVertex->GetXaxis()->SetTitle("z");
474  fHistZVertex->GetYaxis()->SetTitle("counts");
475  fOutput->Add(fHistZVertex);
476 
477  if (fForceBeamType != kpp) {
478  fHistCentrality = new TH1F("fHistCentrality","Event centrality distribution", 200, 0, 100);
479  fHistCentrality->GetXaxis()->SetTitle("Centrality (%)");
480  fHistCentrality->GetYaxis()->SetTitle("counts");
481  fOutput->Add(fHistCentrality);
482 
483  fHistEventPlane = new TH1F("fHistEventPlane","Event plane", 120, -TMath::Pi(), TMath::Pi());
484  fHistEventPlane->GetXaxis()->SetTitle("event plane");
485  fHistEventPlane->GetYaxis()->SetTitle("counts");
486  fOutput->Add(fHistEventPlane);
487  }
488 
490  fHistEventRejection = new TH1F("fHistEventRejection","Reasons to reject event",20,0,20);
491 #if ROOT_VERSION_CODE < ROOT_VERSION(6,4,2)
492  fHistEventRejection->SetBit(TH1::kCanRebin);
493 #else
494  fHistEventRejection->SetCanExtend(TH1::kAllAxes);
495 #endif
496  fHistEventRejection->GetXaxis()->SetBinLabel(1,"PhysSel");
497  fHistEventRejection->GetXaxis()->SetBinLabel(2,"trigger");
498  fHistEventRejection->GetXaxis()->SetBinLabel(3,"trigTypeSel");
499  fHistEventRejection->GetXaxis()->SetBinLabel(4,"Cent");
500  fHistEventRejection->GetXaxis()->SetBinLabel(5,"vertex contr.");
501  fHistEventRejection->GetXaxis()->SetBinLabel(6,"Vz");
502  fHistEventRejection->GetXaxis()->SetBinLabel(7,"VzSPD");
503  fHistEventRejection->GetXaxis()->SetBinLabel(8,"trackInEmcal");
504  fHistEventRejection->GetXaxis()->SetBinLabel(9,"minNTrack");
505  fHistEventRejection->GetXaxis()->SetBinLabel(10,"VtxSel2013pA");
506  fHistEventRejection->GetXaxis()->SetBinLabel(11,"PileUp");
507  fHistEventRejection->GetXaxis()->SetBinLabel(12,"EvtPlane");
508  fHistEventRejection->GetXaxis()->SetBinLabel(13,"SelPtHardBin");
509  fHistEventRejection->GetXaxis()->SetBinLabel(14,"Bkg evt");
510  fHistEventRejection->GetXaxis()->SetBinLabel(15,"RecycleEmbeddedEvent");
511  fHistEventRejection->GetYaxis()->SetTitle("counts");
513  }
514 
515  fHistTriggerClasses = new TH1F("fHistTriggerClasses","fHistTriggerClasses",3,0,3);
516 #if ROOT_VERSION_CODE < ROOT_VERSION(6,4,2)
517  fHistTriggerClasses->SetBit(TH1::kCanRebin);
518 #else
519  fHistTriggerClasses->SetCanExtend(TH1::kAllAxes);
520 #endif
522 
524  fHistTriggerClassesCorr = new TH1F("fHistTriggerClassesCorr","fHistTriggerClassesCorr",3,0,3);
525 #if ROOT_VERSION_CODE < ROOT_VERSION(6,4,2)
526  fHistTriggerClassesCorr->SetBit(TH1::kCanRebin);
527 #else
528  fHistTriggerClassesCorr->SetCanExtend(TH1::kAllAxes);
529 #endif
531  }
532 
533  fHistEventCount = new TH1F("fHistEventCount","fHistEventCount",2,0,2);
534  fHistEventCount->GetXaxis()->SetBinLabel(1,"Accepted");
535  fHistEventCount->GetXaxis()->SetBinLabel(2,"Rejected");
536  fHistEventCount->GetYaxis()->SetTitle("counts");
537  fOutput->Add(fHistEventCount);
538 
539  PostData(1, fOutput);
540 }
541 
543 {
544  if (fIsPythia || fIsHerwig) {
545  // Protection: In case the pt-hard bin handling is not initialized we fall back to the
546  // global pt-hard bin (usually 0) in order to aviod mismatch between histograms before
547  // and after selection
551  fHistPtHard->Fill(fPtHard);
552  if(fPtHardInitialized){
553  fHistPtHardCorr->Fill(fPtHardBin, fPtHard);
555  fHistPtHardBinCorr->Fill(fPtHardBin, fPtHardBinGlobal);
556  }
557  }
558 
559 
560  fHistZVertex->Fill(fVertex[2]);
561 
562  if (fForceBeamType != kpp) {
563  fHistCentrality->Fill(fCent);
564  fHistEventPlane->Fill(fEPV0);
565  }
566 
567  std::unique_ptr<TObjArray> triggerClasses(InputEvent()->GetFiredTriggerClasses().Tokenize(" "));
568  TObjString* triggerClass(nullptr);
569  for(auto trg : *triggerClasses){
570  triggerClass = static_cast<TObjString*>(trg);
571  fHistTriggerClasses->Fill(triggerClass->GetString(), 1);
572  }
573 
575  // downscale-corrected number of events are calculated based on the min. bias reference
576  // Formula: N_corr = N_MB * d_Trg/d_{Min_Bias}
577  if(InputEvent()->GetFiredTriggerClasses().Contains(fMinBiasRefTrigger)){
579  Double_t downscaleref = downscalefactors->GetDownscaleFactorForTriggerClass(fMinBiasRefTrigger);
580  for(auto t : downscalefactors->GetTriggerClasses()){
581  Double_t downscaletrg = downscalefactors->GetDownscaleFactorForTriggerClass(t);
582  fHistTriggerClassesCorr->Fill(t, downscaletrg/downscaleref);
583  }
584  }
585  }
586 
587  return kTRUE;
588 }
589 
591 {
592  // Recycle embedded events which do not pass the internal event selection in the embedding helper
594  auto embeddingHelper = AliAnalysisTaskEmcalEmbeddingHelper::GetInstance();
595  if (embeddingHelper && embeddingHelper->EmbeddedEventUsed() == false) {
596  if (fGeneralHistograms) {
597  fHistEventRejection->Fill("RecycleEmbeddedEvent",1);
598  fHistEventCount->Fill("Rejected",1);
599  }
600  return;
601  }
602  }
603 
604  if (!fLocalInitialized){
605  ExecOnce();
606  UserExecOnce();
607 
608  // Initialize event cuts here: This prevents a segfault
609  // in case a user task overwrites the function UserCreateOutputObjects
611  fAliEventCuts = new AliEventCuts(true);
612  // Do not perform trigger selection in the AliEvent cuts but let the task do this before
613  fAliEventCuts->OverrideAutomaticTriggerSelection(AliVEvent::kAny, true);
614  if(fOutput) fOutput->Add(fAliEventCuts);
615  }
616  }
617 
618  if (!fLocalInitialized)
619  return;
620 
621  if(fFileChanged){
622  FileChanged();
623  fFileChanged = kFALSE;
624  }
625 
626  if (!RetrieveEventObjects())
627  return;
628 
629  if(InputEvent()->GetRunNumber() != fRunNumber){
630  fRunNumber = InputEvent()->GetRunNumber();
633  }
634 
635  // Apply fallback for pythia cross section if needed
637  AliDebugStream(1) << "Fallback to cross section from pythia header required" << std::endl;
638  /*
639  // Get the pthard bin
640  Float_t pthard = fPythiaHeader->GetPtHard();
641  int pthardbin = 0;
642  if(fPtHardBinning.GetSize()){
643  for(int ib = 0; ib < fNPtHardBins; ib++){
644  if(pthard >= static_cast<Float_t>(fPtHardBinning[ib]) && pthard < static_cast<Float_t>(fPtHardBinning[ib+1])) {
645  pthardbin = ib;
646  break;
647  }
648  }
649  }
650  */
651  fHistXsection->Fill(fPtHardBinGlobal, fPythiaHeader->GetXsection());
652  }
653 
654  if (IsEventSelected()) {
655  if (fGeneralHistograms) fHistEventCount->Fill("Accepted",1);
656  }
657  else {
658  if (fGeneralHistograms) fHistEventCount->Fill("Rejected",1);
659  return;
660  }
661 
663  if (!FillGeneralHistograms())
664  return;
665  }
666 
667  if (!Run())
668  return;
669 
670  if (fCreateHisto) {
671  if (!FillHistograms())
672  return;
673  }
674 
675  if (fCreateHisto && fOutput) {
676  // information for this iteration of the UserExec in the container
677  PostData(1, fOutput);
678  }
679 }
680 
682 {
683  AliWarning("AliAnalysisTaskEmcal::AcceptCluster method is deprecated. Please use GetCusterContainer(c)->AcceptCluster(clus).");
684 
685  if (!clus) return kFALSE;
686 
688  if (!cont) {
689  AliError(Form("%s:Container %d not found",GetName(),c));
690  return 0;
691  }
692  UInt_t rejectionReason = 0;
693  return cont->AcceptCluster(clus, rejectionReason);
694 }
695 
696 Bool_t AliAnalysisTaskEmcal::AcceptTrack(AliVParticle *track, Int_t c) const
697 {
698  AliWarning("AliAnalysisTaskEmcal::AcceptTrack method is deprecated. Please use GetParticleContainer(c)->AcceptParticle(clus).");
699 
700  if (!track) return kFALSE;
701 
703  if (!cont) {
704  AliError(Form("%s:Container %d not found",GetName(),c));
705  return 0;
706  }
707 
708  UInt_t rejectionReason = 0;
709  return cont->AcceptParticle(track, rejectionReason);
710 }
711 
712 Bool_t AliAnalysisTaskEmcal::PythiaInfoFromFile(const char* currFile, Float_t &fXsec, Float_t &fTrials, Int_t &pthard)
713 {
714  TString file(currFile);
715  fXsec = 0;
716  fTrials = 1;
717 
718  // Determine archive type
719  TString archivetype;
720  std::unique_ptr<TObjArray> walk(file.Tokenize("/"));
721  for(auto t : *walk){
722  TString &tok = static_cast<TObjString *>(t)->String();
723  if(tok.Contains(".zip")){
724  archivetype = tok;
725  Int_t pos = archivetype.Index(".zip");
726  archivetype.Replace(pos, archivetype.Length() - pos, "");
727  }
728  }
729  if(archivetype.Length()){
730  AliDebugStream(1) << "Auto-detected archive type " << archivetype << std::endl;
731  Ssiz_t pos1 = file.Index(archivetype,archivetype.Length(),0,TString::kExact);
732  Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
733  Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
734  file.Replace(pos+1,pos2-pos1,"");
735  } else {
736  // not an archive take the basename....
737  file.ReplaceAll(gSystem->BaseName(file.Data()),"");
738  }
739  AliDebugStream(1) << "File name: " << file << std::endl;
740 
741  // Build virtual file name
742  // Support for train tests
743  TString virtualFileName;
744  if(file.Contains("__alice")){
745  TString tmp(file);
746  Int_t pos = tmp.Index("__alice");
747  tmp.Replace(0, pos, "");
748  tmp.ReplaceAll("__", "/");
749  // cut out tag for archive and root file
750  // this needs a determin
751  std::unique_ptr<TObjArray> toks(tmp.Tokenize("/"));
752  TString tag = "_" + archivetype;
753  for(auto t : *toks){
754  TString &path = static_cast<TObjString *>(t)->String();
755  if(path.Contains(tag)){
756  Int_t posTag = path.Index(tag);
757  path.Replace(posTag, path.Length() - posTag, "");
758  }
759  virtualFileName += "/" + path;
760  }
761  } else {
762  virtualFileName = file;
763  }
764 
765  AliDebugStream(1) << "Physical file name " << file << ", virtual file name " << virtualFileName << std::endl;
766 
767  // Get the pt hard bin
768  TString strPthard(virtualFileName);
769 
770  /*
771  // Dead code - to be removed after testing phase
772  // Procedure will fail for everything else than the expected path name
773  strPthard.Remove(strPthard.Last('/'));
774  strPthard.Remove(strPthard.Last('/'));
775  if (strPthard.Contains("AOD")) strPthard.Remove(strPthard.Last('/'));
776  strPthard.Remove(0,strPthard.Last('/')+1);
777  if (strPthard.IsDec()) pthard = strPthard.Atoi();
778  else
779  AliWarningStream() << "Could not extract file number from path " << strPthard << std::endl;
780  */
781 
782  // New implementation : pattern matching
783  // Reason: Implementation valid only for old productions (new productions swap run number and pt-hard bin)
784  // Idea: Don't use the position in the string but the match different informations
785  // + Year clearly 2000+
786  // + Run number can be match to the one in the event
787  // + If we know it is not year or run number, it must be the pt-hard bin if we start from the beginning
788  // The procedure is only valid for the current implementations and unable to detect non-pt-hard bins
789  // It will also fail in case of arbitrary file names
790 
791  bool binfound = false;
792  std::unique_ptr<TObjArray> tokens(strPthard.Tokenize("/"));
793  for(auto t : *tokens) {
794  TString &tok = static_cast<TObjString *>(t)->String();
795  if(tok.IsDec()){
796  Int_t number = tok.Atoi();
797  if(number > 2000 && number < 3000){
798  // Year
799  continue;
800  } else if(number == fInputHandler->GetEvent()->GetRunNumber()){
801  // Run number
802  continue;
803  } else {
804  if(!binfound){
805  // the first number that is not one of the two must be the pt-hard bin
806  binfound = true;
807  pthard = number;
808  break;
809  }
810  }
811  }
812  }
813  if(!binfound) {
814  AliErrorStream() << "Could not extract file number from path " << strPthard << std::endl;
815  } else {
816  AliInfoStream() << "Auto-detecting pt-hard bin " << pthard << std::endl;
817  }
818 
819  AliInfoStream() << "File: " << file << std::endl;
820 
821  // 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
822  std::unique_ptr<TFile> fxsec(TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")));
823 
824  if (!fxsec) {
825  // next trial fetch the histgram file
826  fxsec = std::unique_ptr<TFile>(TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root")));
827  if (!fxsec){
828  AliErrorStream() << "Failed reading cross section from file " << file << std::endl;
829  fUseXsecFromHeader = true;
830  return kFALSE; // not a severe condition but inciate that we have no information
831  }
832  else {
833  // find the tlist we want to be independtent of the name so use the Tkey
834  TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
835  if (!key) return kFALSE;
836  TList *list = dynamic_cast<TList*>(key->ReadObj());
837  if (!list) return kFALSE;
838  TProfile *xSecHist = static_cast<TProfile*>(list->FindObject("h1Xsec"));
839  // check for failure
840  if(!xSecHist->GetEntries()) {
841  // No cross seciton information available - fall back to raw
842  AliErrorStream() << "No cross section information available in file " << fxsec->GetName() <<" - fall back to cross section in PYTHIA header" << std::endl;
843  fUseXsecFromHeader = true;
844  } else {
845  // Cross section histogram filled - take it from there
846  fXsec = xSecHist->GetBinContent(1);
847  if(!fXsec) AliErrorStream() << GetName() << ": Cross section 0 for file " << file << std::endl;
848  fUseXsecFromHeader = false;
849  }
850  fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
851  }
852  } else { // no tree pyxsec.root
853  TTree *xtree = (TTree*)fxsec->Get("Xsection");
854  if (!xtree) return kFALSE;
855  UInt_t ntrials = 0;
856  Double_t xsection = 0;
857  xtree->SetBranchAddress("xsection",&xsection);
858  xtree->SetBranchAddress("ntrials",&ntrials);
859  xtree->GetEntry(0);
860  fTrials = ntrials;
861  fXsec = xsection;
862  }
863  return kTRUE;
864 }
865 
867  fPtHardInitialized = kFALSE;
868  fFileChanged = kTRUE;
869  return kTRUE;
870 }
871 
874  return kTRUE;
875 
876  // Debugging:
877  AliInfoStream() << "FileChanged called for run " << InputEvent()->GetRunNumber() << std::endl;
878 
879  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
880  if (!tree) {
881  AliErrorStream() << GetName() << " - FileChanged: No current tree!" << std::endl;
882  return kFALSE;
883  }
884 
885  Float_t xsection = 0;
886  Float_t trials = 0;
887  Int_t pthardbin = -1;
888 
889  TFile *curfile = tree->GetCurrentFile();
890  if (!curfile) {
891  AliErrorStream() << GetName() << " - FileChanged: No current file!" << std::endl;
892  return kFALSE;
893  }
894 
895  TChain *chain = dynamic_cast<TChain*>(tree);
896  if (chain) tree = chain->GetTree();
897 
898  Int_t nevents = tree->GetEntriesFast();
899 
900  fUseXsecFromHeader = false;
901  PythiaInfoFromFile(curfile->GetName(), xsection, trials, pthardbin);
903  // Use the bin obtained from the path name
904  fPtHardBinGlobal = pthardbin;
905  fPtHardInitialized = kTRUE;
906  } else {
907  // Put everything in the first bin
908  fPtHardBinGlobal = 0;
909  pthardbin = 0;
910  }
911 
912  if ((pthardbin < 0) || (pthardbin > fNPtHardBins-1)){
913  AliErrorStream() << GetName() << ": Invalid global pt-hard bin " << pthardbin << " detected" << std::endl;
914  pthardbin = 0;
915  }
916  fHistTrials->Fill(pthardbin, trials);
917  if(!fUseXsecFromHeader){
918  AliDebugStream(1) << "Using cross section from file pyxsec.root" << std::endl;
919  fHistXsection->Fill(pthardbin, xsection);
920  }
921  fHistEvents->Fill(pthardbin, nevents);
922 
923  return kTRUE;
924 }
925 
927 {
928  if (!fPythiaInfoName.IsNull() && !fPythiaInfo) {
929  fPythiaInfo = dynamic_cast<AliEmcalPythiaInfo*>(event->FindListObject(fPythiaInfoName));
930  if (!fPythiaInfo) {
931  AliError(Form("%s: Could not retrieve parton infos! %s!", GetName(), fPythiaInfoName.Data()));
932  return;
933  }
934  }
935 }
936 
938 {
939  if (!InputEvent()) {
940  AliError(Form("%s: Could not retrieve event! Returning!", GetName()));
941  return;
942  }
943 
944  LoadPythiaInfo(InputEvent());
945 
946  if (fNeedEmcalGeom) {
947  fGeom = AliEMCALGeometry::GetInstanceFromRunNumber(InputEvent()->GetRunNumber());
948  if (!fGeom) {
949  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()));
950  return;
951  }
952  }
953 
954  if (fEventPlaneVsEmcal >= 0) {
955  if (fGeom) {
956  Double_t ep = (fGeom->GetArm1PhiMax() + fGeom->GetArm1PhiMin()) / 2 * TMath::DegToRad() + fEventPlaneVsEmcal - TMath::Pi();
957  fMinEventPlane = ep - TMath::Pi() / 4;
958  fMaxEventPlane = ep + TMath::Pi() / 4;
959  }
960  else {
961  AliWarning("Could not set event plane limits because EMCal geometry was not loaded!");
962  }
963  }
964 
965  //Load all requested track branches - each container knows name already
966  for (Int_t i =0; i<fParticleCollArray.GetEntriesFast(); i++) {
967  AliParticleContainer *cont = static_cast<AliParticleContainer*>(fParticleCollArray.At(i));
968  cont->SetArray(InputEvent());
969  }
970 
971  if (fParticleCollArray.GetEntriesFast()>0) {
973  if (!fTracks) {
974  AliError(Form("%s: Could not retrieve first track branch!", GetName()));
975  return;
976  }
977  }
978 
979  //Load all requested cluster branches - each container knows name already
980  for (Int_t i =0; i<fClusterCollArray.GetEntriesFast(); i++) {
981  AliClusterContainer *cont = static_cast<AliClusterContainer*>(fClusterCollArray.At(i));
982  cont->SetArray(InputEvent());
983  }
984 
985  if (fClusterCollArray.GetEntriesFast()>0) {
987  if (!fCaloClusters) {
988  AliError(Form("%s: Could not retrieve first cluster branch!", GetName()));
989  return;
990  }
991  }
992 
993  if (!fCaloCellsName.IsNull() && !fCaloCells) {
994  fCaloCells = dynamic_cast<AliVCaloCells*>(InputEvent()->FindListObject(fCaloCellsName));
995  if (!fCaloCells) {
996  AliError(Form("%s: Could not retrieve cells %s!", GetName(), fCaloCellsName.Data()));
997  return;
998  }
999  }
1000 
1001  if (!fCaloTriggersName.IsNull() && !fCaloTriggers) {
1002  fCaloTriggers = dynamic_cast<AliVCaloTrigger*>(InputEvent()->FindListObject(fCaloTriggersName));
1003  if (!fCaloTriggers) {
1004  AliError(Form("%s: Could not retrieve calo triggers %s!", GetName(), fCaloTriggersName.Data()));
1005  return;
1006  }
1007  }
1008 
1009  if (!fCaloTriggerPatchInfoName.IsNull() && !fTriggerPatchInfo) {
1010  fTriggerPatchInfo = GetArrayFromEvent(fCaloTriggerPatchInfoName.Data(),"AliEMCALTriggerPatchInfo");
1011  if (!fTriggerPatchInfo) {
1012  AliError(Form("%s: Could not retrieve calo trigger patch info %s!", GetName(), fCaloTriggerPatchInfoName.Data()));
1013  return;
1014  }
1015 
1016  }
1017 
1018  fLocalInitialized = kTRUE;
1019 }
1020 
1022 {
1023  if (fForceBeamType != kNA)
1024  return fForceBeamType;
1025 
1026  AliESDEvent *esd = dynamic_cast<AliESDEvent*>(InputEvent());
1027  if (esd) {
1028  const AliESDRun *run = esd->GetESDRun();
1029  TString beamType = run->GetBeamType();
1030  if (beamType == "p-p")
1031  return kpp;
1032  else if (beamType == "A-A")
1033  return kAA;
1034  else if (beamType == "p-A")
1035  return kpA;
1036  else
1037  return kNA;
1038  } else {
1039  Int_t runNumber = InputEvent()->GetRunNumber();
1040  // All run number ranges taken from the RCT
1041  if ((runNumber >= 136833 && runNumber <= 139517) || // LHC10h
1042  (runNumber >= 167693 && runNumber <= 170593) || // LHC11h
1043  (runNumber >= 244824 && runNumber <= 246994)) { // LHC15o
1044  return kAA;
1045  } else if ((runNumber >= 188356 && runNumber <= 188366) || // LHC12g
1046  (runNumber >= 195164 && runNumber <= 197388) || // LHC13b-f
1047  (runNumber >= 265015 && runNumber <= 267166)) { // LHC16q-t
1048  return kpA;
1049  } else {
1050  return kpp;
1051  }
1052  }
1053 }
1054 
1056 {
1057  if (!fTriggerPatchInfo)
1058  return 0;
1059 
1060  //number of patches in event
1061  Int_t nPatch = fTriggerPatchInfo->GetEntries();
1062 
1063  //loop over patches to define trigger type of event
1064  Int_t nG1 = 0;
1065  Int_t nG2 = 0;
1066  Int_t nJ1 = 0;
1067  Int_t nJ2 = 0;
1068  Int_t nL0 = 0;
1069  AliEMCALTriggerPatchInfo *patch;
1070  for (Int_t iPatch = 0; iPatch < nPatch; iPatch++) {
1071  patch = (AliEMCALTriggerPatchInfo*)fTriggerPatchInfo->At( iPatch );
1072  if (patch->IsGammaHigh()) nG1++;
1073  if (patch->IsGammaLow()) nG2++;
1074  if (patch->IsJetHigh()) nJ1++;
1075  if (patch->IsJetLow()) nJ2++;
1076  if (patch->IsLevel0()) nL0++;
1077  }
1078 
1079  AliDebug(2, "Patch summary: ");
1080  AliDebug(2, Form("Number of patches: %d", nPatch));
1081  AliDebug(2, Form("Jet: low[%d], high[%d]" ,nJ2, nJ1));
1082  AliDebug(2, Form("Gamma: low[%d], high[%d]" ,nG2, nG1));
1083 
1084  ULong_t triggers(0);
1085  if (nL0>0) SETBIT(triggers, kL0);
1086  if (nG1>0) SETBIT(triggers, kG1);
1087  if (nG2>0) SETBIT(triggers, kG2);
1088  if (nJ1>0) SETBIT(triggers, kJ1);
1089  if (nJ2>0) SETBIT(triggers, kJ2);
1090  return triggers;
1091 }
1092 
1094 {
1095  //
1096  if(trigger==kND) {
1097  AliWarning(Form("%s: Requesting undefined trigger type!", GetName()));
1098  return kFALSE;
1099  }
1100  //MV: removing this logic which as far as I can see doesn't make any sense
1101  // if(trigger & kND){
1102  // return fTriggers == 0;
1103  // }
1104  return TESTBIT(fTriggers, trigger);
1105 }
1106 
1109  if(!IsTriggerSelected()) return false;
1110  if(!CheckMCOutliers()) return false;
1111  return fAliEventCuts->AcceptEvent(fInputEvent);
1112 }
1113 
1115 {
1116  AliDebugStream(1) << "Using default event selection" << std::endl;
1117  if (fOffTrigger != AliVEvent::kAny) {
1118  UInt_t res = 0;
1119  const AliESDEvent *eev = dynamic_cast<const AliESDEvent*>(InputEvent());
1120  if (eev) {
1121  res = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
1122  } else {
1123  const AliAODEvent *aev = dynamic_cast<const AliAODEvent*>(InputEvent());
1124  if (aev) {
1125  res = ((AliVAODHeader*)aev->GetHeader())->GetOfflineTrigger();
1126  }
1127  }
1128  if ((res & fOffTrigger) == 0) {
1129  if (fGeneralHistograms) fHistEventRejection->Fill("PhysSel",1);
1130  return kFALSE;
1131  }
1132  }
1133 
1134  if(!IsTriggerSelected()) {
1135  if (fGeneralHistograms) fHistEventRejection->Fill("trigger",1);
1136  return kFALSE;
1137  }
1138 
1139  if (fTriggerTypeSel != kND) {
1141  if (fGeneralHistograms) fHistEventRejection->Fill("trigTypeSel",1);
1142  return kFALSE;
1143  }
1144  }
1145 
1146  if ((fMinCent != -999) && (fMaxCent != -999)) {
1147  if (fCent<fMinCent || fCent>fMaxCent) {
1148  if (fGeneralHistograms) fHistEventRejection->Fill("Cent",1);
1149  return kFALSE;
1150  }
1151  }
1152 
1153  if (fUseAliAnaUtils) {
1154  if (!fAliAnalysisUtils)
1155  fAliAnalysisUtils = new AliAnalysisUtils();
1156  fAliAnalysisUtils->SetMinVtxContr(fMinVertexContrib);
1157  fAliAnalysisUtils->SetMaxVtxZ(999);
1158  if(fMinVz<-998.) fMinVz = -10.;
1159  if(fMaxVz>998.) fMaxVz = 10.;
1160 
1161  if (!fAliAnalysisUtils->IsVertexSelected2013pA(InputEvent())) {
1162  if (fGeneralHistograms) fHistEventRejection->Fill("VtxSel2013pA",1);
1163  return kFALSE;
1164  }
1165 
1166  if (fRejectPileup && fAliAnalysisUtils->IsPileUpEvent(InputEvent())) {
1167  if (fGeneralHistograms) fHistEventRejection->Fill("PileUp",1);
1168  return kFALSE;
1169  }
1170 
1171  if(fTklVsClusSPDCut && fAliAnalysisUtils->IsSPDClusterVsTrackletBG(InputEvent())) {
1172  if (fGeneralHistograms) fHistEventRejection->Fill("Bkg evt",1);
1173  return kFALSE;
1174  }
1175  }
1176 
1177  if ((fMinVz > -998.) && (fMaxVz < 998.)) {
1178  if (fNVertCont == 0 ) {
1179  if (fGeneralHistograms) fHistEventRejection->Fill("vertex contr.",1);
1180  return kFALSE;
1181  }
1182  Double_t vz = fVertex[2];
1183  if (vz < fMinVz || vz > fMaxVz) {
1184  if (fGeneralHistograms) fHistEventRejection->Fill("Vz",1);
1185  return kFALSE;
1186  }
1187 
1188  if (fNVertSPDCont > 0 && fZvertexDiff < 999) {
1189  Double_t vzSPD = fVertexSPD[2];
1190  Double_t dvertex = TMath::Abs(vz-vzSPD);
1191  //if difference larger than fZvertexDiff
1192  if (dvertex > fZvertexDiff) {
1193  if (fGeneralHistograms) fHistEventRejection->Fill("VzSPD",1);
1194  return kFALSE;
1195  }
1196  }
1197  }
1198 
1199  if (fMinPtTrackInEmcal > 0 && fGeom) {
1200  Bool_t trackInEmcalOk = kFALSE;
1201  Int_t ntracks = GetNParticles(0);
1202  for (Int_t i = 0; i < ntracks; i++) {
1203  AliVParticle *track = GetAcceptParticleFromArray(i,0);
1204  if (!track)
1205  continue;
1206 
1207  Double_t phiMin = fGeom->GetArm1PhiMin() * TMath::DegToRad();
1208  Double_t phiMax = fGeom->GetArm1PhiMax() * TMath::DegToRad();
1209  Int_t runNumber = InputEvent()->GetRunNumber();
1210  if (runNumber>=177295 && runNumber<=197470) { //small SM masked in 2012 and 2013
1211  phiMin = 1.4;
1212  phiMax = TMath::Pi();
1213  }
1214 
1215  if (track->Eta() < fGeom->GetArm1EtaMin() || track->Eta() > fGeom->GetArm1EtaMax() || track->Phi() < phiMin || track->Phi() > phiMax)
1216  continue;
1217  if (track->Pt() > fMinPtTrackInEmcal) {
1218  trackInEmcalOk = kTRUE;
1219  break;
1220  }
1221  }
1222  if (!trackInEmcalOk) {
1223  if (fGeneralHistograms) fHistEventRejection->Fill("trackInEmcal",1);
1224  return kFALSE;
1225  }
1226  }
1227 
1228  if (fMinNTrack > 0) {
1229  Int_t nTracksAcc = 0;
1230  Int_t ntracks = GetNParticles(0);
1231  for (Int_t i = 0; i < ntracks; i++) {
1232  AliVParticle *track = GetAcceptParticleFromArray(i,0);
1233  if (!track)
1234  continue;
1235  if (track->Pt() > fTrackPtCut) {
1236  nTracksAcc++;
1237  if (nTracksAcc>=fMinNTrack)
1238  break;
1239  }
1240  }
1241  if (nTracksAcc<fMinNTrack) {
1242  if (fGeneralHistograms) fHistEventRejection->Fill("minNTrack",1);
1243  return kFALSE;
1244  }
1245  }
1246 
1247  if (!(fEPV0 > fMinEventPlane && fEPV0 <= fMaxEventPlane) &&
1248  !(fEPV0 + TMath::Pi() > fMinEventPlane && fEPV0 + TMath::Pi() <= fMaxEventPlane) &&
1249  !(fEPV0 - TMath::Pi() > fMinEventPlane && fEPV0 - TMath::Pi() <= fMaxEventPlane))
1250  {
1251  if (fGeneralHistograms) fHistEventRejection->Fill("EvtPlane",1);
1252  return kFALSE;
1253  }
1254 
1255  if (fSelectPtHardBin != -999 && fSelectPtHardBin != fPtHardBin) {
1256  if (fGeneralHistograms) fHistEventRejection->Fill("SelPtHardBin",1);
1257  return kFALSE;
1258  }
1259 
1260  // Reject filter for MC data
1261  if (!CheckMCOutliers()) return kFALSE;
1262 
1263  return kTRUE;
1264 }
1265 
1267  // Default implementation of trigger selection
1268  // same as previously (code moved from IsEventSelected
1269  // to trigger selection). Users should re-implement
1270  // this function in case they have certain needs, in
1271  // particular for EMCAL triggers
1272  AliDebugStream(1) << "Using default trigger selection" << std::endl;
1273  if (!fTrigClass.IsNull()) {
1274  TString fired = InputEvent()->GetFiredTriggerClasses();
1275  if (!fired.Contains("-B-")) return kFALSE;
1276 
1277  std::unique_ptr<TObjArray> arr(fTrigClass.Tokenize("|"));
1278  if (!arr) return kFALSE;
1279  Bool_t match = false;
1280  for (Int_t i=0;i<arr->GetEntriesFast();++i) {
1281  TObject *obj = arr->At(i);
1282  if (!obj) continue;
1283 
1284  //Check if requested trigger was fired
1285  TString objStr = obj->GetName();
1287  (objStr.Contains("J1") || objStr.Contains("J2") || objStr.Contains("G1") || objStr.Contains("G2"))) {
1288  // This is relevant for EMCal triggers with 2 thresholds
1289  // If the kOverlapWithLowThreshold was requested than the overlap between the two triggers goes with the lower threshold trigger
1290  TString trigType1 = "J1";
1291  TString trigType2 = "J2";
1292  if(objStr.Contains("G")) {
1293  trigType1 = "G1";
1294  trigType2 = "G2";
1295  }
1296  if(objStr.Contains(trigType2) && fired.Contains(trigType2.Data())) { //requesting low threshold + overlap
1297  match = 1;
1298  break;
1299  } else if(objStr.Contains(trigType1) && fired.Contains(trigType1.Data()) && !fired.Contains(trigType2.Data())) { //high threshold only
1300  match = 1;
1301  break;
1302  }
1303  }
1304  else {
1305  // If this is not an EMCal trigger, or no particular treatment of EMCal triggers was requested,
1306  // simply check that the trigger was fired
1307  if (fired.Contains(obj->GetName())) {
1308  match = 1;
1309  break;
1310  }
1311  }
1312  }
1313  if (!match) return kFALSE;
1314  }
1315  return kTRUE;
1316 }
1317 
1319 {
1320  if (!fPythiaHeader || !fMCRejectFilter) return kTRUE;
1321 
1322  // Condition 1: Pythia jet / pT-hard > factor
1323  if (fPtHardAndJetPtFactor > 0.) {
1324  AliTLorentzVector jet;
1325 
1326  Int_t nTriggerJets = fPythiaHeader->NTriggerJets();
1327 
1328  AliDebug(1,Form("Njets: %d, pT Hard %f",nTriggerJets, fPtHard));
1329 
1330  Float_t tmpjet[]={0,0,0,0};
1331  for (Int_t ijet = 0; ijet< nTriggerJets; ijet++) {
1332  fPythiaHeader->TriggerJet(ijet, tmpjet);
1333 
1334  jet.SetPxPyPzE(tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3]);
1335 
1336  AliDebug(1,Form("jet %d; pycell jet pT %f",ijet, jet.Pt()));
1337 
1338  //Compare jet pT and pt Hard
1339  if (jet.Pt() > fPtHardAndJetPtFactor * fPtHard) {
1340  AliInfo(Form("Reject jet event with : pT Hard %2.2f, pycell jet pT %2.2f, rejection factor %1.1f\n", fPtHard, jet.Pt(), fPtHardAndJetPtFactor));
1341  return kFALSE;
1342  }
1343  }
1344  }
1345  // end condition 1
1346 
1347  // Condition 2 : Reconstructed EMCal cluster pT / pT-hard > factor
1348  if (fPtHardAndClusterPtFactor > 0.) {
1349  AliClusterContainer* mccluscont = GetClusterContainer(0);
1350  if ((Bool_t)mccluscont) {
1351  for (auto cluster : mccluscont->all()) {// Not cuts applied ; use accept for cuts
1352  Float_t ecluster = cluster->E();
1353 
1354  if (ecluster > (fPtHardAndClusterPtFactor * fPtHard)) {
1355  AliInfo(Form("Reject : ecluster %2.2f, calo %d, factor %2.2f, ptHard %f",ecluster,cluster->GetType(),fPtHardAndClusterPtFactor,fPtHard));
1356  return kFALSE;
1357  }
1358  }
1359  }
1360  }
1361  // end condition 2
1362 
1363  // condition 3 : Reconstructed track pT / pT-hard >factor
1364  if (fPtHardAndTrackPtFactor > 0.) {
1365  AliMCParticleContainer* mcpartcont = dynamic_cast<AliMCParticleContainer*>(GetParticleContainer(0));
1366  if ((Bool_t)mcpartcont) {
1367  for (auto mctrack : mcpartcont->all()) {// Not cuts applied ; use accept for cuts
1368  Float_t trackpt = mctrack->Pt();
1369  if (trackpt > (fPtHardAndTrackPtFactor * fPtHard) ) {
1370  AliInfo(Form("Reject : track %2.2f, factor %2.2f, ptHard %f", trackpt, fPtHardAndTrackPtFactor, fPtHard));
1371  return kFALSE;
1372  }
1373  }
1374  }
1375  }
1376  // end condition 3
1377 
1378  return kTRUE;
1379 }
1380 
1381 TClonesArray *AliAnalysisTaskEmcal::GetArrayFromEvent(const char *name, const char *clname)
1382 {
1383  TClonesArray *arr = 0;
1384  TString sname(name);
1385  if (!sname.IsNull()) {
1386  arr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(sname));
1387  if (!arr) {
1388  AliWarning(Form("%s: Could not retrieve array with name %s!", GetName(), name));
1389  return 0;
1390  }
1391  } else {
1392  return 0;
1393  }
1394 
1395  if (!clname)
1396  return arr;
1397 
1398  TString objname(arr->GetClass()->GetName());
1399  TClass cls(objname);
1400  if (!cls.InheritsFrom(clname)) {
1401  AliWarning(Form("%s: Objects of type %s in %s are not inherited from %s!",
1402  GetName(), cls.GetName(), name, clname));
1403  return 0;
1404  }
1405  return arr;
1406 }
1407 
1409 {
1410  fVertex[0] = 0;
1411  fVertex[1] = 0;
1412  fVertex[2] = 0;
1413  fNVertCont = 0;
1414 
1415  fVertexSPD[0] = 0;
1416  fVertexSPD[1] = 0;
1417  fVertexSPD[2] = 0;
1418  fNVertSPDCont = 0;
1419 
1420  if (fGeneratePythiaInfoObject && MCEvent()) {
1421  GeneratePythiaInfoObject(MCEvent());
1422  }
1423 
1424  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
1425  if (vert) {
1426  vert->GetXYZ(fVertex);
1427  fNVertCont = vert->GetNContributors();
1428  }
1429 
1430  const AliVVertex *vertSPD = InputEvent()->GetPrimaryVertexSPD();
1431  if (vertSPD) {
1432  vertSPD->GetXYZ(fVertexSPD);
1433  fNVertSPDCont = vertSPD->GetNContributors();
1434  }
1435 
1436  fBeamType = GetBeamType();
1437  TObject * header = InputEvent()->GetHeader();
1438  if (fBeamType == kAA || fBeamType == kpA ) {
1440  if (header->InheritsFrom("AliNanoAODStorage")){
1441  AliNanoAODHeader *nanoHead = (AliNanoAODHeader*)header;
1442  fCent=nanoHead->GetCentr(fCentEst.Data());
1443  }else{
1444  AliMultSelection *MultSelection = static_cast<AliMultSelection*>(InputEvent()->FindListObject("MultSelection"));
1445  if (MultSelection) {
1446  fCent = MultSelection->GetMultiplicityPercentile(fCentEst.Data());
1447  }
1448  else {
1449  AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
1450  }
1451  }
1452  }
1453  else { // old centrality estimation < 2015
1454  if (header->InheritsFrom("AliNanoAODStorage")){
1455  AliNanoAODHeader *nanoHead = (AliNanoAODHeader*)header;
1456  fCent=nanoHead->GetCentr(fCentEst.Data());
1457  }else{
1458  AliCentrality *aliCent = InputEvent()->GetCentrality();
1459  if (aliCent) {
1460  fCent = aliCent->GetCentralityPercentile(fCentEst.Data());
1461  }
1462  else {
1463  AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
1464  }
1465  }
1466  }
1467 
1468  if (fNcentBins==4) {
1469  if (fCent >= 0 && fCent < 10) fCentBin = 0;
1470  else if (fCent >= 10 && fCent < 30) fCentBin = 1;
1471  else if (fCent >= 30 && fCent < 50) fCentBin = 2;
1472  else if (fCent >= 50 && fCent <= 100) fCentBin = 3;
1473  else {
1474  AliWarning(Form("%s: Negative centrality: %f. Assuming 99", GetName(), fCent));
1475  fCentBin = fNcentBins-1;
1476  }
1477  }
1478  else if (fNcentBins==5) { // for PbPb 2015
1479  if (fCent >= 0 && fCent < 10) fCentBin = 0;
1480  else if (fCent >= 10 && fCent < 30) fCentBin = 1;
1481  else if (fCent >= 30 && fCent < 50) fCentBin = 2;
1482  else if (fCent >= 50 && fCent <= 90) fCentBin = 3;
1483  else if (fCent > 90) {
1484  fCent = 99;
1485  fCentBin = 4;
1486  }
1487  else {
1488  AliWarning(Form("%s: Negative centrality: %f. Assuming 99", GetName(), fCent));
1489  fCentBin = fNcentBins-1;
1490  }
1491  }
1492  else {
1493  Double_t centWidth = (fMaxCent-fMinCent)/(Double_t)fNcentBins;
1494  if(centWidth>0.) {
1495  fCentBin = TMath::FloorNint(fCent/centWidth);
1496  }
1497  else {
1498  fCentBin = 0;
1499  }
1500  if (fCentBin>=fNcentBins) {
1501  AliWarning(Form("%s: fCentBin too large: cent = %f fCentBin = %d. Assuming 99", GetName(),fCent,fCentBin));
1502  fCentBin = fNcentBins-1;
1503  }
1504  }
1505  if (header->InheritsFrom("AliNanoAODStorage")){
1506  AliNanoAODHeader *nanoHead = (AliNanoAODHeader*)header;
1507  fEPV0=nanoHead->GetVar(nanoHead->GetVarIndex("cstEvPlaneV0"));
1508  fEPV0A=nanoHead->GetVar(nanoHead->GetVarIndex("cstEvPlaneV0A"));
1509  fEPV0C=nanoHead->GetVar(nanoHead->GetVarIndex("cstEvPlaneV0C"));
1510  }else{
1511  AliEventplane *aliEP = InputEvent()->GetEventplane();
1512  if (aliEP) {
1513  fEPV0 = aliEP->GetEventplane("V0" ,InputEvent());
1514  fEPV0A = aliEP->GetEventplane("V0A",InputEvent());
1515  fEPV0C = aliEP->GetEventplane("V0C",InputEvent());
1516  } else {
1517  AliWarning(Form("%s: Could not retrieve event plane information!", GetName()));
1518  }
1519  }
1520  }
1521  else {
1522  fCent = 99;
1523  fCentBin = 0;
1524  }
1525 
1526  if (fIsPythia) {
1527  if (MCEvent()) {
1528  fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
1529  if (!fPythiaHeader) {
1530  // Check if AOD
1531  AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
1532 
1533  if (aodMCH) {
1534  for (UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++) {
1535  fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
1536  if (fPythiaHeader) break;
1537  }
1538  }
1539  }
1540  }
1541  }
1542 
1543  if (fPythiaHeader) {
1544  fPtHard = fPythiaHeader->GetPtHard();
1545 
1546  if(fPtHardBinning.GetSize()){
1547  // pt-hard binning defined for the corresponding dataset - automatically determine the bin
1548  for (fPtHardBin = 0; fPtHardBin < fNPtHardBins; fPtHardBin++) {
1549  if (fPtHard >= static_cast<float>(fPtHardBinning[fPtHardBin]) && fPtHard < static_cast<float>(fPtHardBinning[fPtHardBin+1]))
1550  break;
1551  }
1552  } else {
1553  // No pt-hard binning defined for the dataset - leaving the bin to 0
1554  fPtHardBin = 0;
1555  }
1556 
1557  if(fPtHardInitialized){
1558  // do check only in case the global pt-hard bin is initialized
1560  AliErrorStream() << GetName() << ": Mismatch in pt-hard bin determination. Local: " << fPtHardBin << ", Global: " << fPtHardBinGlobal << std::endl;
1561  }
1562  }
1563 
1564  fXsection = fPythiaHeader->GetXsection();
1565  fNTrials = fPythiaHeader->Trials();
1566  }
1567 
1568  if (fIsHerwig) {
1569  if (MCEvent()) {
1570  fHerwigHeader = dynamic_cast<AliGenHerwigEventHeader*>(MCEvent()->GenEventHeader());
1571 
1572  if (!fHerwigHeader) {
1573  // Check if AOD
1574  AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
1575 
1576  if (aodMCH) {
1577  for (UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++) {
1578  fHerwigHeader = dynamic_cast<AliGenHerwigEventHeader*>(aodMCH->GetCocktailHeader(i));
1579  if (fHerwigHeader) break;
1580  }
1581  }
1582  }
1583  }
1584  }
1585 
1586  if (fHerwigHeader) {
1587  fPtHard = fHerwigHeader->GetPtHard();
1588 
1589  if(fPtHardBinning.GetSize()){
1590  // pt-hard binning defined for the corresponding dataset - automatically determine the bin
1591  for (fPtHardBin = 0; fPtHardBin < fNPtHardBins; fPtHardBin++) {
1593  break;
1594  }
1595  } else {
1596  // No pt-hard binning defined for the dataset - leaving the bin to 0
1597  fPtHardBin = 0;
1598  }
1599  fXsection = fHerwigHeader->Weight();
1600  fNTrials = fHerwigHeader->Trials();
1601  }
1602 
1603 
1605 
1606  AliEmcalContainer* cont = 0;
1607 
1608  TIter nextPartColl(&fParticleCollArray);
1609  while ((cont = static_cast<AliEmcalContainer*>(nextPartColl()))){
1610  cont->NextEvent(InputEvent());
1611  }
1612 
1613  TIter nextClusColl(&fClusterCollArray);
1614  while ((cont = static_cast<AliParticleContainer*>(nextClusColl()))){
1615  cont->NextEvent(InputEvent());
1616  }
1617 
1618  return kTRUE;
1619 }
1620 
1622 {
1623  if (TString(n).IsNull()) return 0;
1624 
1626 
1627  fParticleCollArray.Add(cont);
1628 
1629  return cont;
1630 }
1631 
1633 {
1634  if (TString(n).IsNull()) return 0;
1635 
1636  AliTrackContainer* cont = new AliTrackContainer(n);
1637 
1638  fParticleCollArray.Add(cont);
1639 
1640  return cont;
1641 }
1642 
1644 {
1645  if (TString(n).IsNull()) return 0;
1646 
1648 
1649  fParticleCollArray.Add(cont);
1650 
1651  return cont;
1652 }
1653 
1655 {
1656  if (TString(n).IsNull()) return 0;
1657 
1659 
1660  fClusterCollArray.Add(cont);
1661 
1662  return cont;
1663 }
1664 
1666 {
1667  if (i<0 || i>fParticleCollArray.GetEntriesFast()) return 0;
1668  AliParticleContainer *cont = static_cast<AliParticleContainer*>(fParticleCollArray.At(i));
1669  return cont;
1670 }
1671 
1673 {
1674  if (i<0 || i>fClusterCollArray.GetEntriesFast()) return 0;
1675  AliClusterContainer *cont = static_cast<AliClusterContainer*>(fClusterCollArray.At(i));
1676  return cont;
1677 }
1678 
1680 {
1681  AliParticleContainer *cont = static_cast<AliParticleContainer*>(fParticleCollArray.FindObject(name));
1682  return cont;
1683 }
1684 
1686 {
1687  AliClusterContainer *cont = static_cast<AliClusterContainer*>(fClusterCollArray.FindObject(name));
1688  return cont;
1689 }
1690 
1692 {
1694  if (!cont) {
1695  AliError(Form("%s: Particle container %d not found",GetName(),i));
1696  return 0;
1697  }
1698  TString contName = cont->GetArrayName();
1699  return cont->GetArray();
1700 }
1701 
1703 {
1705  if (!cont) {
1706  AliError(Form("%s:Cluster container %d not found",GetName(),i));
1707  return 0;
1708  }
1709  return cont->GetArray();
1710 }
1711 
1713 {
1714 
1716  if (!cont) {
1717  AliError(Form("%s: Particle container %d not found",GetName(),c));
1718  return 0;
1719  }
1720  AliVParticle *vp = cont->GetAcceptParticle(p);
1721 
1722  return vp;
1723 }
1724 
1726 {
1728  if (!cont) {
1729  AliError(Form("%s: Cluster container %d not found",GetName(),c));
1730  return 0;
1731  }
1732  AliVCluster *vc = cont->GetAcceptCluster(cl);
1733 
1734  return vc;
1735 }
1736 
1738 {
1740  if (!cont) {
1741  AliError(Form("%s: Particle container %d not found",GetName(),i));
1742  return 0;
1743  }
1744  return cont->GetNEntries();
1745 }
1746 
1748 {
1750  if (!cont) {
1751  AliError(Form("%s: Cluster container %d not found",GetName(),i));
1752  return 0;
1753  }
1754  return cont->GetNEntries();
1755 }
1756 
1757 AliEMCALTriggerPatchInfo* AliAnalysisTaskEmcal::GetMainTriggerPatch(TriggerCategory trigger, Bool_t doSimpleOffline)
1758 {
1759 
1760  if (!fTriggerPatchInfo) {
1761  AliError(Form("%s: fTriggerPatchInfo not available",GetName()));
1762  return 0;
1763  }
1764 
1765  //number of patches in event
1766  Int_t nPatch = fTriggerPatchInfo->GetEntries();
1767 
1768  //extract main trigger patch(es)
1769  AliEMCALTriggerPatchInfo *patch(NULL), *selected(NULL);
1770  for (Int_t iPatch = 0; iPatch < nPatch; iPatch++) {
1771 
1772  patch = (AliEMCALTriggerPatchInfo*)fTriggerPatchInfo->At( iPatch );
1773  if (patch->IsMainTrigger()) {
1774  if(doSimpleOffline){
1775  if(patch->IsOfflineSimple()){
1776  switch(trigger){
1777  case kTriggerLevel0:
1778  // option not yet implemented in the trigger maker
1779  if(patch->IsLevel0()) selected = patch;
1780  break;
1781  case kTriggerLevel1Jet:
1782  if(patch->IsJetHighSimple() || patch->IsJetLowSimple()){
1783  if(!selected) selected = patch;
1784  else if(patch->GetADCOfflineAmp() > selected->GetADCOfflineAmp()) selected = patch;
1785  }
1786  break;
1787  case kTriggerLevel1Gamma:
1788  if(patch->IsGammaHighSimple() || patch->IsGammaLowSimple()){
1789  if(!selected) selected = patch;
1790  else if(patch->GetADCOfflineAmp() > selected->GetADCOfflineAmp()) selected = patch;
1791  }
1792  break;
1793  default: // Silence compiler warnings
1794  AliError("Untreated case: Main Patch is recalculated; should be in 'else' branch");
1795  };
1796  }
1797  } else { // Not OfflineSimple
1798  switch(trigger){
1799  case kTriggerLevel0:
1800  if(patch->IsLevel0()) selected = patch;
1801  break;
1802  case kTriggerLevel1Jet:
1803  if(patch->IsJetHigh() || patch->IsJetLow()){
1804  if(!selected) selected = patch;
1805  else if (patch->GetADCAmp() > selected->GetADCAmp())
1806  selected = patch;
1807  }
1808  break;
1809  case kTriggerLevel1Gamma:
1810  if(patch->IsGammaHigh() || patch->IsGammaLow()){
1811  if(!selected) selected = patch;
1812  else if (patch->GetADCAmp() > selected->GetADCAmp())
1813  selected = patch;
1814  }
1815  break;
1816  default:
1817  AliError("Untreated case: Main Patch is recalculated; should be in 'else' branch");
1818  };
1819  }
1820  }
1821  else if ((trigger == kTriggerRecalcJet && patch->IsRecalcJet()) ||
1822  (trigger == kTriggerRecalcGamma && patch->IsRecalcGamma())) { // recalculated patches
1823  if (doSimpleOffline && patch->IsOfflineSimple()) {
1824  if(!selected) selected = patch;
1825  else if (patch->GetADCOfflineAmp() > selected->GetADCOfflineAmp()) // this in fact should not be needed, but we have it in teh other branches as well, so keeping it for compleness
1826  selected = patch;
1827  }
1828  else if (!doSimpleOffline && !patch->IsOfflineSimple()) {
1829  if(!selected) selected = patch;
1830  else if (patch->GetADCAmp() > selected->GetADCAmp())
1831  selected = patch;
1832  }
1833  }
1834  }
1835  return selected;
1836 }
1837 
1839 {
1840  if (!(InputEvent()->FindListObject(obj->GetName()))) {
1841  InputEvent()->AddObject(obj);
1842  }
1843  else {
1844  if (!attempt) {
1845  AliFatal(Form("%s: Container with name %s already present. Aborting", GetName(), obj->GetName()));
1846  }
1847  }
1848 }
1849 
1851 {
1852 
1853  if (!fGeom) {
1854  AliWarning(Form("%s - AliAnalysisTaskEmcal::IsTrackInEmcalAcceptance - Geometry is not available!", GetName()));
1855  return kFALSE;
1856  }
1857 
1858  Double_t minPhi = fGeom->GetArm1PhiMin() - edges;
1859  Double_t maxPhi = fGeom->GetArm1PhiMax() + edges;
1860 
1861  if (part->Phi() > minPhi && part->Phi() < maxPhi) {
1862  return kTRUE;
1863  }
1864  else {
1865  return kFALSE;
1866  }
1867 }
1868 
1870 {
1871  axis->SetBinLabel(1, "NullObject");
1872  axis->SetBinLabel(2, "Pt");
1873  axis->SetBinLabel(3, "Acceptance");
1874  axis->SetBinLabel(4, "MCLabel");
1875  axis->SetBinLabel(5, "BitMap");
1876  axis->SetBinLabel(6, "HF cut");
1877  axis->SetBinLabel(7, "Bit6");
1878  axis->SetBinLabel(8, "NotHybridTrack");
1879  axis->SetBinLabel(9, "MCFlag");
1880  axis->SetBinLabel(10, "MCGenerator");
1881  axis->SetBinLabel(11, "ChargeCut");
1882  axis->SetBinLabel(12, "MinDistanceTPCSectorEdge");
1883  axis->SetBinLabel(13, "Bit12");
1884  axis->SetBinLabel(14, "IsEMCal");
1885  axis->SetBinLabel(15, "Time");
1886  axis->SetBinLabel(16, "Energy");
1887  axis->SetBinLabel(17, "ExoticCut");
1888  axis->SetBinLabel(18, "Bit17");
1889  axis->SetBinLabel(19, "Area");
1890  axis->SetBinLabel(20, "AreaEmc");
1891  axis->SetBinLabel(21, "ZLeadingCh");
1892  axis->SetBinLabel(22, "ZLeadingEmc");
1893  axis->SetBinLabel(23, "NEF");
1894  axis->SetBinLabel(24, "MinLeadPt");
1895  axis->SetBinLabel(25, "MaxTrackPt");
1896  axis->SetBinLabel(26, "MaxClusterPt");
1897  axis->SetBinLabel(27, "Flavour");
1898  axis->SetBinLabel(28, "TagStatus");
1899  axis->SetBinLabel(29, "MinNConstituents");
1900  axis->SetBinLabel(30, "Bit29");
1901  axis->SetBinLabel(31, "Bit30");
1902  axis->SetBinLabel(32, "Bit31");
1903 }
1904 
1905 Double_t AliAnalysisTaskEmcal::GetParallelFraction(AliVParticle* part1, AliVParticle* part2)
1906 {
1907  TVector3 vect1(part1->Px(), part1->Py(), part1->Pz());
1908  TVector3 vect2(part2->Px(), part2->Py(), part2->Pz());
1909  Double_t z = (vect1 * vect2) / (vect2 * vect2);
1910  return z;
1911 }
1912 
1913 Double_t AliAnalysisTaskEmcal::GetParallelFraction(const TVector3& vect1, AliVParticle* part2)
1914 {
1915  TVector3 vect2(part2->Px(), part2->Py(), part2->Pz());
1916  Double_t z = (vect1 * vect2) / (vect2 * vect2);
1917  return z;
1918 }
1919 
1920 void AliAnalysisTaskEmcal::GetEtaPhiDiff(const AliVTrack *t, const AliVCluster *v, Double_t &phidiff, Double_t &etadiff)
1921 {
1922  phidiff = 999;
1923  etadiff = 999;
1924 
1925  if (!t||!v) return;
1926 
1927  Double_t veta = t->GetTrackEtaOnEMCal();
1928  Double_t vphi = t->GetTrackPhiOnEMCal();
1929 
1930  Float_t pos[3] = {0};
1931  v->GetPosition(pos);
1932  TVector3 cpos(pos);
1933  Double_t ceta = cpos.Eta();
1934  Double_t cphi = cpos.Phi();
1935  etadiff=veta-ceta;
1936  phidiff=TVector2::Phi_mpi_pi(vphi-cphi);
1937 }
1938 
1939 Byte_t AliAnalysisTaskEmcal::GetTrackType(const AliVTrack *t)
1940 {
1941  Byte_t ret = 0;
1942  if (t->TestBit(BIT(22)) && !t->TestBit(BIT(23)))
1943  ret = 1;
1944  else if (!t->TestBit(BIT(22)) && t->TestBit(BIT(23)))
1945  ret = 2;
1946  else if (t->TestBit(BIT(22)) && t->TestBit(BIT(23)))
1947  ret = 3;
1948  return ret;
1949 }
1950 
1951 Byte_t AliAnalysisTaskEmcal::GetTrackType(const AliAODTrack *aodTrack, UInt_t filterBit1, UInt_t filterBit2)
1952 {
1953 
1954  Int_t res = 0;
1955 
1956  if (aodTrack->TestFilterBit(filterBit1)) {
1957  res = 0;
1958  }
1959  else if (aodTrack->TestFilterBit(filterBit2)) {
1960  if ((aodTrack->GetStatus()&AliVTrack::kITSrefit)!=0) {
1961  res = 1;
1962  }
1963  else {
1964  res = 2;
1965  }
1966  }
1967  else {
1968  res = 3;
1969  }
1970 
1971  return res;
1972 }
1973 
1975 {
1976  if (!fPythiaInfo) {
1978  }
1979 
1980  AliStack* stack = mcEvent->Stack();
1981 
1982  const Int_t nprim = stack->GetNprimary();
1983  // reject if partons are missing from stack for some reason
1984  if (nprim < 8) return;
1985 
1986  TParticle *part6 = stack->Particle(6);
1987  TParticle *part7 = stack->Particle(7);
1988 
1989  fPythiaInfo->SetPartonFlag6(TMath::Abs(part6->GetPdgCode()));
1990  fPythiaInfo->SetParton6(part6->Pt(), part6->Eta(), part6->Phi(), part6->GetMass());
1991 
1992  fPythiaInfo->SetPartonFlag7(TMath::Abs(part7->GetPdgCode()));
1993  fPythiaInfo->SetParton7(part7->Pt(), part7->Eta(), part7->Phi(), part7->GetMass());
1994 
1995  AliGenPythiaEventHeader *pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(mcEvent->GenEventHeader());
1996  if(pythiaGenHeader){
1997  Float_t ptWeight=pythiaGenHeader->EventWeight();
1998  fPythiaInfo->SetPythiaEventWeight(ptWeight);}
1999 }
2000 
2002 {
2003  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
2004  if (!mgr) {
2005  ::Error("AddAODHandler", "No analysis manager to connect to.");
2006  return NULL;
2007  }
2008 
2009  AliAODInputHandler* aodHandler = new AliAODInputHandler();
2010 
2011  AliVEventHandler *inputHandler=mgr->GetInputEventHandler();
2012  if (inputHandler && (inputHandler->IsA() == AliMultiInputEventHandler::Class())) {
2013  AliMultiInputEventHandler *multiInputHandler=(AliMultiInputEventHandler*)inputHandler;
2014  multiInputHandler->AddInputEventHandler(aodHandler);
2015  }
2016  else {
2017  if (!inputHandler) {
2018  mgr->SetInputEventHandler(aodHandler);
2019  }
2020  else {
2021  ::Error("AddAODHandler", "inputHandler is NOT null. AOD handler was NOT added !!!");
2022  return NULL;
2023  }
2024  }
2025 
2026  return aodHandler;
2027 }
2028 
2030 {
2031  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
2032  if (!mgr) {
2033  ::Error("AddESDHandler", "No analysis manager to connect to.");
2034  return NULL;
2035  }
2036 
2037  AliESDInputHandler *esdHandler = new AliESDInputHandler();
2038 
2039  AliVEventHandler *inputHandler=mgr->GetInputEventHandler();
2040  if (inputHandler && (inputHandler->IsA() == AliMultiInputEventHandler::Class())) {
2041  AliMultiInputEventHandler *multiInputHandler=(AliMultiInputEventHandler*)inputHandler;
2042  multiInputHandler->AddInputEventHandler(esdHandler);
2043  }
2044  else {
2045  if (!inputHandler) {
2046  mgr->SetInputEventHandler(esdHandler);
2047  }
2048  else {
2049  ::Error("AddESDHandler", "inputHandler is NOT null. ESD handler was NOT added !!!");
2050  return NULL;
2051  }
2052  }
2053 
2054  return esdHandler;
2055 }
2056 
virtual Bool_t Run()
Run function. This is the core function of the analysis and contains the user code. Therefore users have to implement this function.
Bool_t fGeneratePythiaInfoObject
Generate Pythia info object.
TObjArray fClusterCollArray
cluster collection array
Int_t fNVertSPDCont
!event SPD vertex number of contributors
void SetParticlePtCut(Double_t cut)
Bool_t fIsPythia
trigger, if it is a PYTHIA production
double Double_t
Definition: External.C:58
void SetParton7(Float_t pt, Float_t eta, Float_t phi, Float_t mass=0)
TH2 * fHistPtHardBinCorr
!Correlation between global and local (per-event) -hard bin
TH1 * fHistTrials
!trials from pyxsec.root
Bool_t fPtHardInitialized
!flag whether the -hard bin was initialized, purely for internal processing
Definition: External.C:236
void SetArray(const AliVEvent *event)
EMCAL Level1 gamma trigger, low threshold.
AliEmcalPythiaInfo * fPythiaInfo
!event parton info
Bool_t AcceptTrack(AliVParticle *track, Int_t c=0) const
ULong_t GetTriggerList()
Get list of selected triggers of the given event.
Int_t fPtHardBinGlobal
!event -hard bin, detected from filename
EMCAL Level1 jet trigger, low threshold.
Bool_t HasTriggerType(TriggerType triggersel)
Check if event has a given trigger type.
Int_t fNTrials
!event trials
UInt_t fOffTrigger
offline trigger for event selection
Double_t fVertexSPD[3]
!event Svertex
Double_t fMinCent
min centrality for event selection
Double_t fTrackPtCut
cut on track pt in event selection
Bool_t fRecycleUnusedEmbeddedEventsMode
Allows the recycling of embedded events which fail internal event selection. See the embedding helper...
Recalculated jet trigger patch; does not need to be above trigger threshold.
Base task in the EMCAL framework.
Bool_t fLocalInitialized
whether or not the task has been already initialized
Bool_t fUsePtHardBinScaling
Use -hard bin scaling in merging.
TH2 * fHistPtHardCorrGlobal
!Correlation between -hard value and global bin
void SetPartonFlag7(Int_t flag7)
Container with name, TClonesArray and cuts for particles.
Bool_t fUseXsecFromHeader
! Use cross section from header instead of pyxsec.root (purely transient)
TSystem * gSystem
void SetTrackPtCut(Double_t cut, Int_t c=0)
Apply cut on the transverse momentum of all tracks in the track container with index c...
static void GetEtaPhiDiff(const AliVTrack *t, const AliVCluster *v, Double_t &phidiff, Double_t &etadiff)
Calculate and difference between a track (t) and a cluster (c).
Double_t fMinBinPt
min pt in histograms
Double_t fEPV0
!event plane V0
Double_t phiMin
static AliEmcalDownscaleFactorsOCDB * Instance()
Int_t fNPtHardBins
Number of -hard bins in the dataset.
Bool_t fGeneralHistograms
whether or not it should fill some general histograms
Declaration of class AliAnalysisTaskEmcalEmbeddingHelper.
Bool_t AcceptCluster(AliVCluster *clus, Int_t c=0) const
Cluster selection.
TCanvas * c
Definition: TestFitELoss.C:172
virtual void UserExecOnce()
Task initializations handled in user tasks.
Int_t fCentBin
!event centrality bin
virtual Bool_t IsTriggerSelected()
Selection of a hardware trigger.
TH1 * fHistEventsAfterSel
!total number of events per pt hard bin after selection
Double_t GetDownscaleFactorForTriggerClass(const TString &trigger) const
const AliMCParticleIterableContainer all() const
Int_t fMinVertexContrib
Min. number of vertex contributors.
Float_t fPtHardAndClusterPtFactor
Factor between ptHard and cluster pT to reject/accept event.
Double_t fMinPtTrackInEmcal
min pt track in emcal
AliStack * stack
TH1 * fHistEventPlane
!event plane distribution
TH1 * fHistEvents
!total number of events per pt hard bin
void SetClusPtCut(Double_t cut, Int_t c=0)
Apply cut on for all clusters in container with index c.
AliClusterContainer * AddClusterContainer(const char *n)
Create new cluster container and attach it to the task.
Double_t fEPV0C
!event plane V0C
void SetParton6(Float_t pt, Float_t eta, Float_t phi, Float_t mass=0)
TH1 * fHistCentrality
!event centrality distribution
Container for particles within the EMCAL framework.
Bool_t fIsEmbedded
trigger, embedded signal
TObjArray fParticleCollArray
particle/track collection array
BeamType
Switch for the beam type.
void SetTrackEtaLimits(Double_t min, Double_t max, Int_t c=0)
Apply cut on the pseudorapidity of the all tracks in the track container with index c...
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
TProfile * fHistXsectionAfterSel
!x section from pythia header
TriggerType
Switch for EMCAL trigger types.
Bool_t CheckMCOutliers()
Filter the mc tails in pt-hard distributions.
EMCalTriggerMode_t fEMCalTriggerMode
EMCal trigger selection mode.
virtual Bool_t FillHistograms()
Function filling histograms.
Int_t GetNParticles(Int_t i=0) const
Get number of particles in container attached to this task with index i.
TClonesArray * fCaloClusters
!clusters
Bool_t fUseNewCentralityEstimation
Use new centrality estimation (for 2015 data)
Bool_t IsTrackInEmcalAcceptance(AliVParticle *part, Double_t edges=0.9) const
Determines if a track is inside the EMCal acceptance.
int Int_t
Definition: External.C:63
TH1 * fHistTriggerClasses
!number of events in each trigger class
Bool_t fCountDownscaleCorrectedEvents
Count event number corrected for downscaling.
unsigned int UInt_t
Definition: External.C:33
float Float_t
Definition: External.C:68
Double_t fMaxVz
max vertex for event selection
void GeneratePythiaInfoObject(AliMCEvent *mcEvent)
Copy some information about the Pythia event in a PythaInfo object.
AliEMCALGeometry * fGeom
!emcal geometry
The overlap between low and high threshold trigger is assigned to the lower threshold only...
kRecalculated gamma trigger patch; does not need to be above trigger threshold
TString fCaloTriggerPatchInfoName
trigger patch info array name
TString fCaloTriggersName
name of calo triggers collection
TH1 * fHistTriggerClassesCorr
!corrected number of events in each trigger class
Bool_t fIsHerwig
trigger, if it is a HERWIG production
AliGenPythiaEventHeader * fPythiaHeader
!event Pythia header
void SetTrackPhiLimits(Double_t min, Double_t max, Int_t c=0)
Apply cut on azimuthal angle of the all tracks in the track container with index c...
Double_t phiMax
AliParticleContainer * AddParticleContainer(const char *n)
Create new particle container and attach it to the task.
AliAnalysisUtils * fAliAnalysisUtils
!vertex selection (optional)
BeamType fForceBeamType
forced beam type
void SetUseScaling(Bool_t val)
Definition: AliEmcalList.h:31
Int_t fNcentBins
how many centrality bins
AliClusterContainer * GetClusterContainer(Int_t i=0) const
Get cluster container attached to this task.
TriggerType fTriggerTypeSel
trigger type to select based on trigger patches
virtual Bool_t AcceptCluster(Int_t i, UInt_t &rejectionReason) const
virtual Bool_t FillGeneralHistograms()
Filling general histograms.
TString fTrigClass
trigger class name for event selection
Float_t fPtHardAndJetPtFactor
Factor between ptHard and jet pT to reject/accept event.
AliVCluster * GetAcceptCluster(Int_t i) const
TClonesArray * GetParticleArray(Int_t i=0) const
Get TClonesArray with particles.
BeamType GetBeamType() const
Get beam type.
Double_t fMinVz
min vertex for event selection
virtual AliVParticle * GetAcceptParticle(Int_t i=-1) const
BeamType fBeamType
!event beam type
Float_t fPtHardAndTrackPtFactor
Factor between ptHard and track pT to reject/accept event.
Double_t fCent
!event centrality
Double_t fMinEventPlane
minimum event plane value
TString fCaloCellsName
name of calo cell collection
Int_t fMCLabelShift
if MC label > fMCLabelShift, MC label -= fMCLabelShift
Int_t GetNClusters(Int_t i=0) const
Get number of clusters in the cluster container attached to this task with index i.
Int_t fNVertCont
!event vertex number of contributors
Bool_t fMCRejectFilter
enable the filtering of events by tail rejection
Bool_t fGetPtHardBinFromName
Obtain pt-hard bin from file path.
unsigned long ULong_t
Definition: External.C:38
Double_t fZvertexDiff
upper limit for distance between primary and SPD vertex
virtual Bool_t AcceptParticle(const AliVParticle *vp, UInt_t &rejectionReason) const
AliEMCALTriggerPatchInfo * GetMainTriggerPatch(TriggerCategory triggersel=kTriggerLevel1Jet, Bool_t doSimpleOffline=kFALSE)
Get main trigger match.
EMCAL Level0 trigger.
EMCAL Level1 jet trigger, high threshold.
Int_t fSelectPtHardBin
select one pt hard bin for analysis
AliMCParticleContainer * AddMCParticleContainer(const char *n)
Create new container for MC particles and attach it to the task.
static Double_t GetParallelFraction(AliVParticle *part1, AliVParticle *part2)
Calculates the fraction of momentum z of part 1 w.r.t. part 2 in the direction of part 2...
virtual Bool_t RetrieveEventObjects()
Retrieve common objects from event.
Bool_t fRejectPileup
Reject pilup using function AliAnalysisUtils::IsPileUpEvent()
Bool_t fUseInternalEventSelection
Use builtin event selection of the AliAnalysisTaskEmcal instead of AliEventCuts.
TProfile * fHistXsection
!x section from pyxsec.root
Bool_t PythiaInfoFromFile(const char *currFile, Float_t &fXsec, Float_t &fTrials, Int_t &pthard)
Loading PYTHIA information from external cross section file into the task.
void UserExec(Option_t *option)
Event loop, called for each event.
Int_t fMinMCLabel
minimum MC label value for the tracks/clusters being considered MC particles
void SetPartonFlag6(Int_t flag6)
AliVCaloCells * fCaloCells
!cells
TClonesArray * GetArrayFromEvent(const char *name, const char *clname=0)
Read a TClonesArray from event.
Enhanced TList-derived class that implements correct merging for pt_hard binned production.
Definition: AliEmcalList.h:25
Double_t fEventPlaneVsEmcal
select events which have a certain event plane wrt the emcal
Bool_t IsEventSelectedInternal()
Perform event selection (old method)
virtual Bool_t IsEventSelected()
Performing event selection.
Float_t fPtHard
!event -hard
TH1 * fHistPtHard
! -hard distribution
void SetParticleEtaLimits(Double_t min, Double_t max)
AliEmcalList * fOutput
!output list
Handler for downscale factors for various triggers obtained from the OCDB.
Double_t fMaxBinPt
max pt in histograms
Int_t fPtHardBin
!event -hard bin
virtual void RunChanged(Int_t)
Process tasks relevant when a file with a different run number is processed.
TClonesArray * fTracks
!tracks
TH1 * fHistTrialsAfterSel
!total number of trials per pt hard bin after selection
AliGenHerwigEventHeader * fHerwigHeader
!event Herwig header
void LoadPythiaInfo(AliVEvent *event)
Load parton info.
Bool_t fIsEsd
!whether it&#39;s an ESD analysis
static AliAODInputHandler * AddAODHandler()
Add an AOD handler to the analysis manager.
Double_t fVertex[3]
!event vertex
AliTrackContainer * AddTrackContainer(const char *n)
Create new track container and attach it to the task.
Bool_t fCreateHisto
whether or not create histograms
Bool_t UserNotify()
Notifying the user that the input data file has changed and performing steps needed to be done...
Store some informaion about a Pythia eventThis class is used to store some information about a Pythia...
TFile * file
TList with histograms for a given trigger.
TH1 * fHistEventRejection
!book keep reasons for rejecting event
Int_t nevents[nsamples]
TClonesArray * fTriggerPatchInfo
!trigger patch info array
TClonesArray * GetClusterArray(Int_t i=0) const
Get TClonesArray with EMCAL clusters.
Bool_t FileChanged()
Steps to be executed when a few file is loaded into the input handler.
Double_t fEPV0A
!event plane V0A
AliEventCuts * fAliEventCuts
!Event cuts (run2 defaults)
virtual void ExecOnce()
Perform steps needed to initialize the analysis.
TString fCentEst
name of V0 centrality estimator
void SetArray(const AliVEvent *event)
TString fPythiaInfoName
name of pythia info object
Int_t GetRunNumber(TString)
Definition: PlotMuonQA.C:2235
Declaration of class AliEmcalPythiaInfo.
const char Option_t
Definition: External.C:48
void SetClusPtCut(Double_t cut)
Int_t fRunNumber
!run number (triggering RunChanged())
EMCAL Level1 gamma trigger, high threshold.
void AddObjectToEvent(TObject *obj, Bool_t attempt=kFALSE)
Add object to event.
AliVCaloTrigger * fCaloTriggers
!calo triggers
void SetRejectionReasonLabels(TAxis *axis)
void UserCreateOutputObjects()
Main initialization function on the worker.
Bool_t fFileChanged
! Signal triggered when the file has changed
TH1 * fHistZVertex
!z vertex position
Int_t fMinNTrack
minimum nr of tracks in event with pT>fTrackPtCut
static Byte_t GetTrackType(const AliVTrack *t)
Get track type encoded from bits 20 and 21.
Bool_t fUseAliAnaUtils
used for LHC13* data: z-vtx, Ncontributors, z-vtx resolution cuts
void SetClusTimeCut(Double_t min, Double_t max, Int_t c=0)
Apply cut on cluster time for clusters in container with index c.
bool Bool_t
Definition: External.C:53
ULong_t fTriggers
list of fired triggers
static AliESDInputHandler * AddESDHandler()
Add a ESD handler to the analysis manager.
Double_t fMaxEventPlane
maximum event plane value
void SetPythiaEventWeight(Float_t ptWeight)
Float_t fXsection
!x-section from pythia header
TH1 * fHistEventCount
!incoming and selected events
Double_t fMaxCent
max centrality for event selection
void SetClusTimeCut(Double_t min, Double_t max)
TriggerCategory
Online trigger categories.
void SetParticlePhiLimits(Double_t min, Double_t max)
AliVParticle * GetAcceptParticleFromArray(Int_t p, Int_t c=0) const
Get particle p if accepted from container with index c If particle not accepted return 0...
Container structure for EMCAL clusters.
TString fMinBiasRefTrigger
Name of the minmum bias reference trigger, used in the calculation of downscale-corrected event numbe...
Container for MC-true particles within the EMCAL framework.
Bool_t fNeedEmcalGeom
whether or not the task needs the emcal geometry
AliVCluster * GetAcceptClusterFromArray(Int_t cl, Int_t c=0) const
Get cluster cl if accepted from container c If particle not accepted return 0.
Int_t fNbins
no. of pt bins
Bool_t fTklVsClusSPDCut
Apply tracklet-vs-cluster SPD cut to reject background events in pp.
AliAnalysisTaskEmcal()
Default constructor.
TH2 * fHistPtHardCorr
!Correlation between -hard value and bin
TArrayI fPtHardBinning
-hard binning
TList * OpenFile(const char *fname)
Definition: DrawAnaELoss.C:65
virtual ~AliAnalysisTaskEmcal()
Destructor.
static const AliAnalysisTaskEmcalEmbeddingHelper * GetInstance()
static Double_t fgkEMCalDCalPhiDivide
phi value used to distinguish between DCal and EMCal