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