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