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