AliPhysics  5364b50 (5364b50)
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
539  fHistTrialsAfterSel->Fill(fPtHardInitialized ? fPtHardBin : fPtHardBinGlobal, fNTrials);
540  fHistXsectionAfterSel->Fill(fPtHardInitialized ? fPtHardBin : fPtHardBinGlobal, fXsection);
541  fHistPtHard->Fill(fPtHard);
542  if(fPtHardInitialized){
544  fHistPtHardCorrGlobal->Fill(fPtHardBinGlobal, 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  // Get the pthard bin
620  Float_t pthard = fPythiaHeader->GetPtHard();
621  int pthardbin = 0;
622  if(fPtHardBinning.GetSize()){
623  for(int ib = 0; ib < fNPtHardBins; ib++){
624  if(pthard >= fPtHardBinning[ib] && pthard < fPtHardBinning[ib+1]) {
625  pthardbin = ib;
626  break;
627  }
628  }
629  }
630  fHistXsection->Fill(pthardbin, fPythiaHeader->GetXsection());
631  }
632 
633  if (IsEventSelected()) {
634  if (fGeneralHistograms) fHistEventCount->Fill("Accepted",1);
635  }
636  else {
637  if (fGeneralHistograms) fHistEventCount->Fill("Rejected",1);
638  return;
639  }
640 
642  if (!FillGeneralHistograms())
643  return;
644  }
645 
646  if (!Run())
647  return;
648 
649  if (fCreateHisto) {
650  if (!FillHistograms())
651  return;
652  }
653 
654  if (fCreateHisto && fOutput) {
655  // information for this iteration of the UserExec in the container
656  PostData(1, fOutput);
657  }
658 }
659 
661 {
662  AliWarning("AliAnalysisTaskEmcal::AcceptCluster method is deprecated. Please use GetCusterContainer(c)->AcceptCluster(clus).");
663 
664  if (!clus) return kFALSE;
665 
667  if (!cont) {
668  AliError(Form("%s:Container %d not found",GetName(),c));
669  return 0;
670  }
671  UInt_t rejectionReason = 0;
672  return cont->AcceptCluster(clus, rejectionReason);
673 }
674 
675 Bool_t AliAnalysisTaskEmcal::AcceptTrack(AliVParticle *track, Int_t c) const
676 {
677  AliWarning("AliAnalysisTaskEmcal::AcceptTrack method is deprecated. Please use GetParticleContainer(c)->AcceptParticle(clus).");
678 
679  if (!track) return kFALSE;
680 
682  if (!cont) {
683  AliError(Form("%s:Container %d not found",GetName(),c));
684  return 0;
685  }
686 
687  UInt_t rejectionReason = 0;
688  return cont->AcceptParticle(track, rejectionReason);
689 }
690 
691 Bool_t AliAnalysisTaskEmcal::PythiaInfoFromFile(const char* currFile, Float_t &fXsec, Float_t &fTrials, Int_t &pthard)
692 {
693  TString file(currFile);
694  fXsec = 0;
695  fTrials = 1;
696 
697  // Determine archive type
698  TString archivetype;
699  std::unique_ptr<TObjArray> walk(file.Tokenize("/"));
700  for(auto t : *walk){
701  TString &tok = static_cast<TObjString *>(t)->String();
702  if(tok.Contains(".zip")){
703  archivetype = tok;
704  Int_t pos = archivetype.Index(".zip");
705  archivetype.Replace(pos, archivetype.Length() - pos, "");
706  }
707  }
708  if(archivetype.Length()){
709  AliDebugStream(1) << "Auto-detected archive type " << archivetype << std::endl;
710  Ssiz_t pos1 = file.Index(archivetype,archivetype.Length(),0,TString::kExact);
711  Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
712  Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
713  file.Replace(pos+1,pos2-pos1,"");
714  } else {
715  // not an archive take the basename....
716  file.ReplaceAll(gSystem->BaseName(file.Data()),"");
717  }
718  AliDebugStream(1) << "File name: " << file << std::endl;
719 
720  // Build virtual file name
721  // Support for train tests
722  TString virtualFileName;
723  if(file.Contains("__alice")){
724  TString tmp(file);
725  Int_t pos = tmp.Index("__alice");
726  tmp.Replace(0, pos, "");
727  tmp.ReplaceAll("__", "/");
728  // cut out tag for archive and root file
729  // this needs a determin
730  std::unique_ptr<TObjArray> toks(tmp.Tokenize("/"));
731  TString tag = "_" + archivetype;
732  for(auto t : *toks){
733  TString &path = static_cast<TObjString *>(t)->String();
734  if(path.Contains(tag)){
735  Int_t posTag = path.Index(tag);
736  path.Replace(posTag, path.Length() - posTag, "");
737  }
738  virtualFileName += "/" + path;
739  }
740  } else {
741  virtualFileName = file;
742  }
743 
744  AliDebugStream(1) << "Physical file name " << file << ", virtual file name " << virtualFileName << std::endl;
745 
746  // Get the pt hard bin
747  TString strPthard(virtualFileName);
748 
749  /*
750  // Dead code - to be removed after testing phase
751  // Procedure will fail for everything else than the expected path name
752  strPthard.Remove(strPthard.Last('/'));
753  strPthard.Remove(strPthard.Last('/'));
754  if (strPthard.Contains("AOD")) strPthard.Remove(strPthard.Last('/'));
755  strPthard.Remove(0,strPthard.Last('/')+1);
756  if (strPthard.IsDec()) pthard = strPthard.Atoi();
757  else
758  AliWarningStream() << "Could not extract file number from path " << strPthard << std::endl;
759  */
760 
761  // New implementation : pattern matching
762  // Reason: Implementation valid only for old productions (new productions swap run number and pt-hard bin)
763  // Idea: Don't use the position in the string but the match different informations
764  // + Year clearly 2000+
765  // + Run number can be match to the one in the event
766  // + If we know it is not year or run number, it must be the pt-hard bin if we start from the beginning
767  // The procedure is only valid for the current implementations and unable to detect non-pt-hard bins
768  // It will also fail in case of arbitrary file names
769 
770  bool binfound = false;
771  std::unique_ptr<TObjArray> tokens(strPthard.Tokenize("/"));
772  for(auto t : *tokens) {
773  TString &tok = static_cast<TObjString *>(t)->String();
774  if(tok.IsDec()){
775  Int_t number = tok.Atoi();
776  if(number > 2000 && number < 3000){
777  // Year
778  continue;
779  } else if(number == fInputHandler->GetEvent()->GetRunNumber()){
780  // Run number
781  continue;
782  } else {
783  if(!binfound){
784  // the first number that is not one of the two must be the pt-hard bin
785  binfound = true;
786  pthard = number;
787  break;
788  }
789  }
790  }
791  }
792  if(!binfound) {
793  AliErrorStream() << "Could not extract file number from path " << strPthard << std::endl;
794  } else {
795  AliInfoStream() << "Auto-detecting pt-hard bin " << pthard << std::endl;
796  }
797 
798  AliInfoStream() << "File: " << file << std::endl;
799 
800  // 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
801  std::unique_ptr<TFile> fxsec(TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")));
802 
803  if (!fxsec) {
804  // next trial fetch the histgram file
805  fxsec = std::unique_ptr<TFile>(TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root")));
806  if (!fxsec){
807  AliErrorStream() << "Failed reading cross section from file " << file << std::endl;
808  fUseXsecFromHeader = true;
809  return kFALSE; // not a severe condition but inciate that we have no information
810  }
811  else {
812  // find the tlist we want to be independtent of the name so use the Tkey
813  TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
814  if (!key) return kFALSE;
815  TList *list = dynamic_cast<TList*>(key->ReadObj());
816  if (!list) return kFALSE;
817  TProfile *xSecHist = static_cast<TProfile*>(list->FindObject("h1Xsec"));
818  // check for failure
819  if(!xSecHist->GetEntries()) {
820  // No cross seciton information available - fall back to raw
821  AliErrorStream() << "No cross section information available in file " << fxsec->GetName() <<" - fall back to cross section in PYTHIA header" << std::endl;
822  fUseXsecFromHeader = true;
823  } else {
824  // Cross section histogram filled - take it from there
825  fXsec = xSecHist->GetBinContent(1);
826  if(!fXsec) AliErrorStream() << GetName() << ": Cross section 0 for file " << file << std::endl;
827  fUseXsecFromHeader = false;
828  }
829  fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
830  }
831  } else { // no tree pyxsec.root
832  TTree *xtree = (TTree*)fxsec->Get("Xsection");
833  if (!xtree) return kFALSE;
834  UInt_t ntrials = 0;
835  Double_t xsection = 0;
836  xtree->SetBranchAddress("xsection",&xsection);
837  xtree->SetBranchAddress("ntrials",&ntrials);
838  xtree->GetEntry(0);
839  fTrials = ntrials;
840  fXsec = xsection;
841  }
842  return kTRUE;
843 }
844 
846  fPtHardInitialized = kFALSE;
847  fFileChanged = kTRUE;
848  return kTRUE;
849 }
850 
853  return kTRUE;
854 
855  // Debugging:
856  AliInfoStream() << "FileChanged called for run " << InputEvent()->GetRunNumber() << std::endl;
857 
858  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
859  if (!tree) {
860  AliErrorStream() << GetName() << " - FileChanged: No current tree!" << std::endl;
861  return kFALSE;
862  }
863 
864  Float_t xsection = 0;
865  Float_t trials = 0;
866  Int_t pthardbin = -1;
867 
868  TFile *curfile = tree->GetCurrentFile();
869  if (!curfile) {
870  AliErrorStream() << GetName() << " - FileChanged: No current file!" << std::endl;
871  return kFALSE;
872  }
873 
874  TChain *chain = dynamic_cast<TChain*>(tree);
875  if (chain) tree = chain->GetTree();
876 
877  Int_t nevents = tree->GetEntriesFast();
878 
879  fUseXsecFromHeader = false;
880  PythiaInfoFromFile(curfile->GetName(), xsection, trials, pthardbin);
882  // Use the bin obtained from the path name
883  fPtHardBinGlobal = pthardbin;
884  fPtHardInitialized = kTRUE;
885  } else {
886  // Put everything in the first bin
887  fPtHardBinGlobal = 0;
888  pthardbin = 0;
889  }
890 
891  if ((pthardbin < 0) || (pthardbin > fNPtHardBins-1)){
892  AliErrorStream() << GetName() << ": Invalid global pt-hard bin " << pthardbin << " detected" << std::endl;
893  pthardbin = 0;
894  }
895  fHistTrials->Fill(pthardbin, trials);
896  if(!fUseXsecFromHeader){
897  AliDebugStream(1) << "Using cross section from file pyxsec.root" << std::endl;
898  fHistXsection->Fill(pthardbin, xsection);
899  }
900  fHistEvents->Fill(pthardbin, nevents);
901 
902  return kTRUE;
903 }
904 
906 {
907  if (!fPythiaInfoName.IsNull() && !fPythiaInfo) {
908  fPythiaInfo = dynamic_cast<AliEmcalPythiaInfo*>(event->FindListObject(fPythiaInfoName));
909  if (!fPythiaInfo) {
910  AliError(Form("%s: Could not retrieve parton infos! %s!", GetName(), fPythiaInfoName.Data()));
911  return;
912  }
913  }
914 }
915 
917 {
918  if (!InputEvent()) {
919  AliError(Form("%s: Could not retrieve event! Returning!", GetName()));
920  return;
921  }
922 
923  LoadPythiaInfo(InputEvent());
924 
925  if (fNeedEmcalGeom) {
926  fGeom = AliEMCALGeometry::GetInstanceFromRunNumber(InputEvent()->GetRunNumber());
927  if (!fGeom) {
928  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()));
929  return;
930  }
931  }
932 
933  if (fEventPlaneVsEmcal >= 0) {
934  if (fGeom) {
935  Double_t ep = (fGeom->GetArm1PhiMax() + fGeom->GetArm1PhiMin()) / 2 * TMath::DegToRad() + fEventPlaneVsEmcal - TMath::Pi();
936  fMinEventPlane = ep - TMath::Pi() / 4;
937  fMaxEventPlane = ep + TMath::Pi() / 4;
938  }
939  else {
940  AliWarning("Could not set event plane limits because EMCal geometry was not loaded!");
941  }
942  }
943 
944  //Load all requested track branches - each container knows name already
945  for (Int_t i =0; i<fParticleCollArray.GetEntriesFast(); i++) {
946  AliParticleContainer *cont = static_cast<AliParticleContainer*>(fParticleCollArray.At(i));
947  cont->SetArray(InputEvent());
948  }
949 
950  if (fParticleCollArray.GetEntriesFast()>0) {
952  if (!fTracks) {
953  AliError(Form("%s: Could not retrieve first track branch!", GetName()));
954  return;
955  }
956  }
957 
958  //Load all requested cluster branches - each container knows name already
959  for (Int_t i =0; i<fClusterCollArray.GetEntriesFast(); i++) {
960  AliClusterContainer *cont = static_cast<AliClusterContainer*>(fClusterCollArray.At(i));
961  cont->SetArray(InputEvent());
962  }
963 
964  if (fClusterCollArray.GetEntriesFast()>0) {
966  if (!fCaloClusters) {
967  AliError(Form("%s: Could not retrieve first cluster branch!", GetName()));
968  return;
969  }
970  }
971 
972  if (!fCaloCellsName.IsNull() && !fCaloCells) {
973  fCaloCells = dynamic_cast<AliVCaloCells*>(InputEvent()->FindListObject(fCaloCellsName));
974  if (!fCaloCells) {
975  AliError(Form("%s: Could not retrieve cells %s!", GetName(), fCaloCellsName.Data()));
976  return;
977  }
978  }
979 
980  if (!fCaloTriggersName.IsNull() && !fCaloTriggers) {
981  fCaloTriggers = dynamic_cast<AliVCaloTrigger*>(InputEvent()->FindListObject(fCaloTriggersName));
982  if (!fCaloTriggers) {
983  AliError(Form("%s: Could not retrieve calo triggers %s!", GetName(), fCaloTriggersName.Data()));
984  return;
985  }
986  }
987 
988  if (!fCaloTriggerPatchInfoName.IsNull() && !fTriggerPatchInfo) {
989  fTriggerPatchInfo = GetArrayFromEvent(fCaloTriggerPatchInfoName.Data(),"AliEMCALTriggerPatchInfo");
990  if (!fTriggerPatchInfo) {
991  AliError(Form("%s: Could not retrieve calo trigger patch info %s!", GetName(), fCaloTriggerPatchInfoName.Data()));
992  return;
993  }
994 
995  }
996 
997  fLocalInitialized = kTRUE;
998 }
999 
1001 {
1002  if (fForceBeamType != kNA)
1003  return fForceBeamType;
1004 
1005  AliESDEvent *esd = dynamic_cast<AliESDEvent*>(InputEvent());
1006  if (esd) {
1007  const AliESDRun *run = esd->GetESDRun();
1008  TString beamType = run->GetBeamType();
1009  if (beamType == "p-p")
1010  return kpp;
1011  else if (beamType == "A-A")
1012  return kAA;
1013  else if (beamType == "p-A")
1014  return kpA;
1015  else
1016  return kNA;
1017  } else {
1018  Int_t runNumber = InputEvent()->GetRunNumber();
1019  // All run number ranges taken from the RCT
1020  if ((runNumber >= 136833 && runNumber <= 139517) || // LHC10h
1021  (runNumber >= 167693 && runNumber <= 170593) || // LHC11h
1022  (runNumber >= 244824 && runNumber <= 246994)) { // LHC15o
1023  return kAA;
1024  } else if ((runNumber >= 188356 && runNumber <= 188366) || // LHC12g
1025  (runNumber >= 195164 && runNumber <= 197388) || // LHC13b-f
1026  (runNumber >= 265015 && runNumber <= 267166)) { // LHC16q-t
1027  return kpA;
1028  } else {
1029  return kpp;
1030  }
1031  }
1032 }
1033 
1035 {
1036  if (!fTriggerPatchInfo)
1037  return 0;
1038 
1039  //number of patches in event
1040  Int_t nPatch = fTriggerPatchInfo->GetEntries();
1041 
1042  //loop over patches to define trigger type of event
1043  Int_t nG1 = 0;
1044  Int_t nG2 = 0;
1045  Int_t nJ1 = 0;
1046  Int_t nJ2 = 0;
1047  Int_t nL0 = 0;
1048  AliEMCALTriggerPatchInfo *patch;
1049  for (Int_t iPatch = 0; iPatch < nPatch; iPatch++) {
1050  patch = (AliEMCALTriggerPatchInfo*)fTriggerPatchInfo->At( iPatch );
1051  if (patch->IsGammaHigh()) nG1++;
1052  if (patch->IsGammaLow()) nG2++;
1053  if (patch->IsJetHigh()) nJ1++;
1054  if (patch->IsJetLow()) nJ2++;
1055  if (patch->IsLevel0()) nL0++;
1056  }
1057 
1058  AliDebug(2, "Patch summary: ");
1059  AliDebug(2, Form("Number of patches: %d", nPatch));
1060  AliDebug(2, Form("Jet: low[%d], high[%d]" ,nJ2, nJ1));
1061  AliDebug(2, Form("Gamma: low[%d], high[%d]" ,nG2, nG1));
1062 
1063  ULong_t triggers(0);
1064  if (nL0>0) SETBIT(triggers, kL0);
1065  if (nG1>0) SETBIT(triggers, kG1);
1066  if (nG2>0) SETBIT(triggers, kG2);
1067  if (nJ1>0) SETBIT(triggers, kJ1);
1068  if (nJ2>0) SETBIT(triggers, kJ2);
1069  return triggers;
1070 }
1071 
1073 {
1074  //
1075  if(trigger==kND) {
1076  AliWarning(Form("%s: Requesting undefined trigger type!", GetName()));
1077  return kFALSE;
1078  }
1079  //MV: removing this logic which as far as I can see doesn't make any sense
1080  // if(trigger & kND){
1081  // return fTriggers == 0;
1082  // }
1083  return TESTBIT(fTriggers, trigger);
1084 }
1085 
1087 {
1088  if (fOffTrigger != AliVEvent::kAny) {
1089  UInt_t res = 0;
1090  const AliESDEvent *eev = dynamic_cast<const AliESDEvent*>(InputEvent());
1091  if (eev) {
1092  res = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
1093  } else {
1094  const AliAODEvent *aev = dynamic_cast<const AliAODEvent*>(InputEvent());
1095  if (aev) {
1096  res = ((AliVAODHeader*)aev->GetHeader())->GetOfflineTrigger();
1097  }
1098  }
1099  if ((res & fOffTrigger) == 0) {
1100  if (fGeneralHistograms) fHistEventRejection->Fill("PhysSel",1);
1101  return kFALSE;
1102  }
1103  }
1104 
1105  if (!fTrigClass.IsNull()) {
1106  TString fired;
1107  const AliESDEvent *eev = dynamic_cast<const AliESDEvent*>(InputEvent());
1108  if (eev) {
1109  fired = eev->GetFiredTriggerClasses();
1110  } else {
1111  const AliAODEvent *aev = dynamic_cast<const AliAODEvent*>(InputEvent());
1112  if (aev) {
1113  fired = aev->GetFiredTriggerClasses();
1114  }
1115  }
1116  if (!fired.Contains("-B-")) {
1117  if (fGeneralHistograms) fHistEventRejection->Fill("trigger",1);
1118  return kFALSE;
1119  }
1120 
1121  std::unique_ptr<TObjArray> arr(fTrigClass.Tokenize("|"));
1122  if (!arr) {
1123  if (fGeneralHistograms) fHistEventRejection->Fill("trigger",1);
1124  return kFALSE;
1125  }
1126  Bool_t match = 0;
1127  for (Int_t i=0;i<arr->GetEntriesFast();++i) {
1128  TObject *obj = arr->At(i);
1129  if (!obj)
1130  continue;
1131 
1132  //Check if requested trigger was fired
1133  TString objStr = obj->GetName();
1135  (objStr.Contains("J1") || objStr.Contains("J2") || objStr.Contains("G1") || objStr.Contains("G2"))) {
1136  // This is relevant for EMCal triggers with 2 thresholds
1137  // If the kOverlapWithLowThreshold was requested than the overlap between the two triggers goes with the lower threshold trigger
1138  TString trigType1 = "J1";
1139  TString trigType2 = "J2";
1140  if(objStr.Contains("G")) {
1141  trigType1 = "G1";
1142  trigType2 = "G2";
1143  }
1144  if(objStr.Contains(trigType2) && fired.Contains(trigType2.Data())) { //requesting low threshold + overlap
1145  match = 1;
1146  break;
1147  }
1148  else if(objStr.Contains(trigType1) && fired.Contains(trigType1.Data()) && !fired.Contains(trigType2.Data())) { //high threshold only
1149  match = 1;
1150  break;
1151  }
1152  }
1153  else {
1154  // If this is not an EMCal trigger, or no particular treatment of EMCal triggers was requested,
1155  // simply check that the trigger was fired
1156  if (fired.Contains(obj->GetName())) {
1157  match = 1;
1158  break;
1159  }
1160  }
1161  }
1162  if (!match) {
1163  if (fGeneralHistograms) fHistEventRejection->Fill("trigger",1);
1164  return kFALSE;
1165  }
1166  }
1167 
1168  if (fTriggerTypeSel != kND) {
1170  if (fGeneralHistograms) fHistEventRejection->Fill("trigTypeSel",1);
1171  return kFALSE;
1172  }
1173  }
1174 
1175  if ((fMinCent != -999) && (fMaxCent != -999)) {
1176  if (fCent<fMinCent || fCent>fMaxCent) {
1177  if (fGeneralHistograms) fHistEventRejection->Fill("Cent",1);
1178  return kFALSE;
1179  }
1180  }
1181 
1182  if (fUseAliAnaUtils) {
1183  if (!fAliAnalysisUtils)
1184  fAliAnalysisUtils = new AliAnalysisUtils();
1185  fAliAnalysisUtils->SetMinVtxContr(2);
1186  fAliAnalysisUtils->SetMaxVtxZ(999);
1187  if(fMinVz<-998.) fMinVz = -10.;
1188  if(fMaxVz>998.) fMaxVz = 10.;
1189 
1190  if (!fAliAnalysisUtils->IsVertexSelected2013pA(InputEvent())) {
1191  if (fGeneralHistograms) fHistEventRejection->Fill("VtxSel2013pA",1);
1192  return kFALSE;
1193  }
1194 
1195  if (fRejectPileup && fAliAnalysisUtils->IsPileUpEvent(InputEvent())) {
1196  if (fGeneralHistograms) fHistEventRejection->Fill("PileUp",1);
1197  return kFALSE;
1198  }
1199 
1200  if(fTklVsClusSPDCut && fAliAnalysisUtils->IsSPDClusterVsTrackletBG(InputEvent())) {
1201  if (fGeneralHistograms) fHistEventRejection->Fill("Bkg evt",1);
1202  return kFALSE;
1203  }
1204  }
1205 
1206  if ((fMinVz > -998.) && (fMaxVz < 998.)) {
1207  if (fNVertCont == 0 ) {
1208  if (fGeneralHistograms) fHistEventRejection->Fill("vertex contr.",1);
1209  return kFALSE;
1210  }
1211  Double_t vz = fVertex[2];
1212  if (vz < fMinVz || vz > fMaxVz) {
1213  if (fGeneralHistograms) fHistEventRejection->Fill("Vz",1);
1214  return kFALSE;
1215  }
1216 
1217  if (fNVertSPDCont > 0 && fZvertexDiff < 999) {
1218  Double_t vzSPD = fVertexSPD[2];
1219  Double_t dvertex = TMath::Abs(vz-vzSPD);
1220  //if difference larger than fZvertexDiff
1221  if (dvertex > fZvertexDiff) {
1222  if (fGeneralHistograms) fHistEventRejection->Fill("VzSPD",1);
1223  return kFALSE;
1224  }
1225  }
1226  }
1227 
1228  if (fMinPtTrackInEmcal > 0 && fGeom) {
1229  Bool_t trackInEmcalOk = kFALSE;
1230  Int_t ntracks = GetNParticles(0);
1231  for (Int_t i = 0; i < ntracks; i++) {
1232  AliVParticle *track = GetAcceptParticleFromArray(i,0);
1233  if (!track)
1234  continue;
1235 
1236  Double_t phiMin = fGeom->GetArm1PhiMin() * TMath::DegToRad();
1237  Double_t phiMax = fGeom->GetArm1PhiMax() * TMath::DegToRad();
1238  Int_t runNumber = InputEvent()->GetRunNumber();
1239  if (runNumber>=177295 && runNumber<=197470) { //small SM masked in 2012 and 2013
1240  phiMin = 1.4;
1241  phiMax = TMath::Pi();
1242  }
1243 
1244  if (track->Eta() < fGeom->GetArm1EtaMin() || track->Eta() > fGeom->GetArm1EtaMax() || track->Phi() < phiMin || track->Phi() > phiMax)
1245  continue;
1246  if (track->Pt() > fMinPtTrackInEmcal) {
1247  trackInEmcalOk = kTRUE;
1248  break;
1249  }
1250  }
1251  if (!trackInEmcalOk) {
1252  if (fGeneralHistograms) fHistEventRejection->Fill("trackInEmcal",1);
1253  return kFALSE;
1254  }
1255  }
1256 
1257  if (fMinNTrack > 0) {
1258  Int_t nTracksAcc = 0;
1259  Int_t ntracks = GetNParticles(0);
1260  for (Int_t i = 0; i < ntracks; i++) {
1261  AliVParticle *track = GetAcceptParticleFromArray(i,0);
1262  if (!track)
1263  continue;
1264  if (track->Pt() > fTrackPtCut) {
1265  nTracksAcc++;
1266  if (nTracksAcc>=fMinNTrack)
1267  break;
1268  }
1269  }
1270  if (nTracksAcc<fMinNTrack) {
1271  if (fGeneralHistograms) fHistEventRejection->Fill("minNTrack",1);
1272  return kFALSE;
1273  }
1274  }
1275 
1276  if (!(fEPV0 > fMinEventPlane && fEPV0 <= fMaxEventPlane) &&
1277  !(fEPV0 + TMath::Pi() > fMinEventPlane && fEPV0 + TMath::Pi() <= fMaxEventPlane) &&
1278  !(fEPV0 - TMath::Pi() > fMinEventPlane && fEPV0 - TMath::Pi() <= fMaxEventPlane))
1279  {
1280  if (fGeneralHistograms) fHistEventRejection->Fill("EvtPlane",1);
1281  return kFALSE;
1282  }
1283 
1284  if (fSelectPtHardBin != -999 && fSelectPtHardBin != fPtHardBin) {
1285  if (fGeneralHistograms) fHistEventRejection->Fill("SelPtHardBin",1);
1286  return kFALSE;
1287  }
1288 
1289  // Reject filter for MC data
1290  if (!CheckMCOutliers()) return kFALSE;
1291 
1292  return kTRUE;
1293 }
1294 
1296 {
1297  if (!fPythiaHeader || !fMCRejectFilter) return kTRUE;
1298 
1299  // Condition 1: Pythia jet / pT-hard > factor
1300  if (fPtHardAndJetPtFactor > 0.) {
1301  AliTLorentzVector jet;
1302 
1303  Int_t nTriggerJets = fPythiaHeader->NTriggerJets();
1304 
1305  AliDebug(1,Form("Njets: %d, pT Hard %f",nTriggerJets, fPtHard));
1306 
1307  Float_t tmpjet[]={0,0,0,0};
1308  for (Int_t ijet = 0; ijet< nTriggerJets; ijet++) {
1309  fPythiaHeader->TriggerJet(ijet, tmpjet);
1310 
1311  jet.SetPxPyPzE(tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3]);
1312 
1313  AliDebug(1,Form("jet %d; pycell jet pT %f",ijet, jet.Pt()));
1314 
1315  //Compare jet pT and pt Hard
1316  if (jet.Pt() > fPtHardAndJetPtFactor * fPtHard) {
1317  AliInfo(Form("Reject jet event with : pT Hard %2.2f, pycell jet pT %2.2f, rejection factor %1.1f\n", fPtHard, jet.Pt(), fPtHardAndJetPtFactor));
1318  return kFALSE;
1319  }
1320  }
1321  }
1322  // end condition 1
1323 
1324  // Condition 2 : Reconstructed EMCal cluster pT / pT-hard > factor
1325  if (fPtHardAndClusterPtFactor > 0.) {
1326  AliClusterContainer* mccluscont = GetClusterContainer(0);
1327  if ((Bool_t)mccluscont) {
1328  for (auto cluster : mccluscont->all()) {// Not cuts applied ; use accept for cuts
1329  Float_t ecluster = cluster->E();
1330 
1331  if (ecluster > (fPtHardAndClusterPtFactor * fPtHard)) {
1332  AliInfo(Form("Reject : ecluster %2.2f, calo %d, factor %2.2f, ptHard %f",ecluster,cluster->GetType(),fPtHardAndClusterPtFactor,fPtHard));
1333  return kFALSE;
1334  }
1335  }
1336  }
1337  }
1338  // end condition 2
1339 
1340  // condition 3 : Reconstructed track pT / pT-hard >factor
1341  if (fPtHardAndTrackPtFactor > 0.) {
1342  AliMCParticleContainer* mcpartcont = dynamic_cast<AliMCParticleContainer*>(GetParticleContainer(0));
1343  if ((Bool_t)mcpartcont) {
1344  for (auto mctrack : mcpartcont->all()) {// Not cuts applied ; use accept for cuts
1345  Float_t trackpt = mctrack->Pt();
1346  if (trackpt > (fPtHardAndTrackPtFactor * fPtHard) ) {
1347  AliInfo(Form("Reject : track %2.2f, factor %2.2f, ptHard %f", trackpt, fPtHardAndTrackPtFactor, fPtHard));
1348  return kFALSE;
1349  }
1350  }
1351  }
1352  }
1353  // end condition 3
1354 
1355  return kTRUE;
1356 }
1357 
1358 TClonesArray *AliAnalysisTaskEmcal::GetArrayFromEvent(const char *name, const char *clname)
1359 {
1360  TClonesArray *arr = 0;
1361  TString sname(name);
1362  if (!sname.IsNull()) {
1363  arr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(sname));
1364  if (!arr) {
1365  AliWarning(Form("%s: Could not retrieve array with name %s!", GetName(), name));
1366  return 0;
1367  }
1368  } else {
1369  return 0;
1370  }
1371 
1372  if (!clname)
1373  return arr;
1374 
1375  TString objname(arr->GetClass()->GetName());
1376  TClass cls(objname);
1377  if (!cls.InheritsFrom(clname)) {
1378  AliWarning(Form("%s: Objects of type %s in %s are not inherited from %s!",
1379  GetName(), cls.GetName(), name, clname));
1380  return 0;
1381  }
1382  return arr;
1383 }
1384 
1386 {
1387  fVertex[0] = 0;
1388  fVertex[1] = 0;
1389  fVertex[2] = 0;
1390  fNVertCont = 0;
1391 
1392  fVertexSPD[0] = 0;
1393  fVertexSPD[1] = 0;
1394  fVertexSPD[2] = 0;
1395  fNVertSPDCont = 0;
1396 
1397  if (fGeneratePythiaInfoObject && MCEvent()) {
1398  GeneratePythiaInfoObject(MCEvent());
1399  }
1400 
1401  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
1402  if (vert) {
1403  vert->GetXYZ(fVertex);
1404  fNVertCont = vert->GetNContributors();
1405  }
1406 
1407  const AliVVertex *vertSPD = InputEvent()->GetPrimaryVertexSPD();
1408  if (vertSPD) {
1409  vertSPD->GetXYZ(fVertexSPD);
1410  fNVertSPDCont = vertSPD->GetNContributors();
1411  }
1412 
1413  fBeamType = GetBeamType();
1414  TObject * header = InputEvent()->GetHeader();
1415  if (fBeamType == kAA || fBeamType == kpA ) {
1417  if (header->InheritsFrom("AliNanoAODStorage")){
1418  AliNanoAODHeader *nanoHead = (AliNanoAODHeader*)header;
1419  fCent=nanoHead->GetCentr(fCentEst.Data());
1420  }else{
1421  AliMultSelection *MultSelection = static_cast<AliMultSelection*>(InputEvent()->FindListObject("MultSelection"));
1422  if (MultSelection) {
1423  fCent = MultSelection->GetMultiplicityPercentile(fCentEst.Data());
1424  }
1425  else {
1426  AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
1427  }
1428  }
1429  }
1430  else { // old centrality estimation < 2015
1431  if (header->InheritsFrom("AliNanoAODStorage")){
1432  AliNanoAODHeader *nanoHead = (AliNanoAODHeader*)header;
1433  fCent=nanoHead->GetCentr(fCentEst.Data());
1434  }else{
1435  AliCentrality *aliCent = InputEvent()->GetCentrality();
1436  if (aliCent) {
1437  fCent = aliCent->GetCentralityPercentile(fCentEst.Data());
1438  }
1439  else {
1440  AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
1441  }
1442  }
1443  }
1444 
1445  if (fNcentBins==4) {
1446  if (fCent >= 0 && fCent < 10) fCentBin = 0;
1447  else if (fCent >= 10 && fCent < 30) fCentBin = 1;
1448  else if (fCent >= 30 && fCent < 50) fCentBin = 2;
1449  else if (fCent >= 50 && fCent <= 100) fCentBin = 3;
1450  else {
1451  AliWarning(Form("%s: Negative centrality: %f. Assuming 99", GetName(), fCent));
1452  fCentBin = fNcentBins-1;
1453  }
1454  }
1455  else if (fNcentBins==5) { // for PbPb 2015
1456  if (fCent >= 0 && fCent < 10) fCentBin = 0;
1457  else if (fCent >= 10 && fCent < 30) fCentBin = 1;
1458  else if (fCent >= 30 && fCent < 50) fCentBin = 2;
1459  else if (fCent >= 50 && fCent <= 90) fCentBin = 3;
1460  else if (fCent > 90) {
1461  fCent = 99;
1462  fCentBin = 4;
1463  }
1464  else {
1465  AliWarning(Form("%s: Negative centrality: %f. Assuming 99", GetName(), fCent));
1466  fCentBin = fNcentBins-1;
1467  }
1468  }
1469  else {
1470  Double_t centWidth = (fMaxCent-fMinCent)/(Double_t)fNcentBins;
1471  if(centWidth>0.) {
1472  fCentBin = TMath::FloorNint(fCent/centWidth);
1473  }
1474  else {
1475  fCentBin = 0;
1476  }
1477  if (fCentBin>=fNcentBins) {
1478  AliWarning(Form("%s: fCentBin too large: cent = %f fCentBin = %d. Assuming 99", GetName(),fCent,fCentBin));
1479  fCentBin = fNcentBins-1;
1480  }
1481  }
1482  if (header->InheritsFrom("AliNanoAODStorage")){
1483  AliNanoAODHeader *nanoHead = (AliNanoAODHeader*)header;
1484  fEPV0=nanoHead->GetVar(nanoHead->GetVarIndex("cstEvPlaneV0"));
1485  fEPV0A=nanoHead->GetVar(nanoHead->GetVarIndex("cstEvPlaneV0A"));
1486  fEPV0C=nanoHead->GetVar(nanoHead->GetVarIndex("cstEvPlaneV0C"));
1487  }else{
1488  AliEventplane *aliEP = InputEvent()->GetEventplane();
1489  if (aliEP) {
1490  fEPV0 = aliEP->GetEventplane("V0" ,InputEvent());
1491  fEPV0A = aliEP->GetEventplane("V0A",InputEvent());
1492  fEPV0C = aliEP->GetEventplane("V0C",InputEvent());
1493  } else {
1494  AliWarning(Form("%s: Could not retrieve event plane information!", GetName()));
1495  }
1496  }
1497  }
1498  else {
1499  fCent = 99;
1500  fCentBin = 0;
1501  }
1502 
1503  if (fIsPythia) {
1504  if (MCEvent()) {
1505  fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
1506  if (!fPythiaHeader) {
1507  // Check if AOD
1508  AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
1509 
1510  if (aodMCH) {
1511  for (UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++) {
1512  fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
1513  if (fPythiaHeader) break;
1514  }
1515  }
1516  }
1517  }
1518  }
1519 
1520  if (fPythiaHeader) {
1521  fPtHard = fPythiaHeader->GetPtHard();
1522 
1523  if(fPtHardBinning.GetSize()){
1524  // pt-hard binning defined for the corresponding dataset - automatically determine the bin
1525  for (fPtHardBin = 0; fPtHardBin < fNPtHardBins; fPtHardBin++) {
1526  if (fPtHard >= static_cast<float>(fPtHardBinning[fPtHardBin]) && fPtHard < static_cast<float>(fPtHardBinning[fPtHardBin+1]))
1527  break;
1528  }
1529  } else {
1530  // No pt-hard binning defined for the dataset - leaving the bin to 0
1531  fPtHardBin = 0;
1532  }
1533 
1534  if(fPtHardInitialized){
1535  // do check only in case the global pt-hard bin is initialized
1537  AliErrorStream() << GetName() << ": Mismatch in pt-hard bin determination. Local: " << fPtHardBin << ", Global: " << fPtHardBinGlobal << std::endl;
1538  }
1539  }
1540 
1541  fXsection = fPythiaHeader->GetXsection();
1542  fNTrials = fPythiaHeader->Trials();
1543  }
1544 
1545  if (fIsHerwig) {
1546  if (MCEvent()) {
1547  fHerwigHeader = dynamic_cast<AliGenHerwigEventHeader*>(MCEvent()->GenEventHeader());
1548 
1549  if (!fHerwigHeader) {
1550  // Check if AOD
1551  AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
1552 
1553  if (aodMCH) {
1554  for (UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++) {
1555  fHerwigHeader = dynamic_cast<AliGenHerwigEventHeader*>(aodMCH->GetCocktailHeader(i));
1556  if (fHerwigHeader) break;
1557  }
1558  }
1559  }
1560  }
1561  }
1562 
1563  if (fHerwigHeader) {
1564  fPtHard = fHerwigHeader->GetPtHard();
1565 
1566  if(fPtHardBinning.GetSize()){
1567  // pt-hard binning defined for the corresponding dataset - automatically determine the bin
1568  for (fPtHardBin = 0; fPtHardBin < fNPtHardBins; fPtHardBin++) {
1570  break;
1571  }
1572  } else {
1573  // No pt-hard binning defined for the dataset - leaving the bin to 0
1574  fPtHardBin = 0;
1575  }
1576  fXsection = fHerwigHeader->Weight();
1577  fNTrials = fHerwigHeader->Trials();
1578  }
1579 
1580 
1582 
1583  AliEmcalContainer* cont = 0;
1584 
1585  TIter nextPartColl(&fParticleCollArray);
1586  while ((cont = static_cast<AliEmcalContainer*>(nextPartColl()))){
1587  cont->NextEvent(InputEvent());
1588  }
1589 
1590  TIter nextClusColl(&fClusterCollArray);
1591  while ((cont = static_cast<AliParticleContainer*>(nextClusColl()))){
1592  cont->NextEvent(InputEvent());
1593  }
1594 
1595  return kTRUE;
1596 }
1597 
1599 {
1600  if (TString(n).IsNull()) return 0;
1601 
1603 
1604  fParticleCollArray.Add(cont);
1605 
1606  return cont;
1607 }
1608 
1610 {
1611  if (TString(n).IsNull()) return 0;
1612 
1613  AliTrackContainer* cont = new AliTrackContainer(n);
1614 
1615  fParticleCollArray.Add(cont);
1616 
1617  return cont;
1618 }
1619 
1621 {
1622  if (TString(n).IsNull()) return 0;
1623 
1625 
1626  fParticleCollArray.Add(cont);
1627 
1628  return cont;
1629 }
1630 
1632 {
1633  if (TString(n).IsNull()) return 0;
1634 
1636 
1637  fClusterCollArray.Add(cont);
1638 
1639  return cont;
1640 }
1641 
1643 {
1644  if (i<0 || i>fParticleCollArray.GetEntriesFast()) return 0;
1645  AliParticleContainer *cont = static_cast<AliParticleContainer*>(fParticleCollArray.At(i));
1646  return cont;
1647 }
1648 
1650 {
1651  if (i<0 || i>fClusterCollArray.GetEntriesFast()) return 0;
1652  AliClusterContainer *cont = static_cast<AliClusterContainer*>(fClusterCollArray.At(i));
1653  return cont;
1654 }
1655 
1657 {
1658  AliParticleContainer *cont = static_cast<AliParticleContainer*>(fParticleCollArray.FindObject(name));
1659  return cont;
1660 }
1661 
1663 {
1664  AliClusterContainer *cont = static_cast<AliClusterContainer*>(fClusterCollArray.FindObject(name));
1665  return cont;
1666 }
1667 
1669 {
1671  if (!cont) {
1672  AliError(Form("%s: Particle container %d not found",GetName(),i));
1673  return 0;
1674  }
1675  TString contName = cont->GetArrayName();
1676  return cont->GetArray();
1677 }
1678 
1680 {
1682  if (!cont) {
1683  AliError(Form("%s:Cluster container %d not found",GetName(),i));
1684  return 0;
1685  }
1686  return cont->GetArray();
1687 }
1688 
1690 {
1691 
1693  if (!cont) {
1694  AliError(Form("%s: Particle container %d not found",GetName(),c));
1695  return 0;
1696  }
1697  AliVParticle *vp = cont->GetAcceptParticle(p);
1698 
1699  return vp;
1700 }
1701 
1703 {
1705  if (!cont) {
1706  AliError(Form("%s: Cluster container %d not found",GetName(),c));
1707  return 0;
1708  }
1709  AliVCluster *vc = cont->GetAcceptCluster(cl);
1710 
1711  return vc;
1712 }
1713 
1715 {
1717  if (!cont) {
1718  AliError(Form("%s: Particle container %d not found",GetName(),i));
1719  return 0;
1720  }
1721  return cont->GetNEntries();
1722 }
1723 
1725 {
1727  if (!cont) {
1728  AliError(Form("%s: Cluster container %d not found",GetName(),i));
1729  return 0;
1730  }
1731  return cont->GetNEntries();
1732 }
1733 
1734 AliEMCALTriggerPatchInfo* AliAnalysisTaskEmcal::GetMainTriggerPatch(TriggerCategory trigger, Bool_t doSimpleOffline)
1735 {
1736 
1737  if (!fTriggerPatchInfo) {
1738  AliError(Form("%s: fTriggerPatchInfo not available",GetName()));
1739  return 0;
1740  }
1741 
1742  //number of patches in event
1743  Int_t nPatch = fTriggerPatchInfo->GetEntries();
1744 
1745  //extract main trigger patch(es)
1746  AliEMCALTriggerPatchInfo *patch(NULL), *selected(NULL);
1747  for (Int_t iPatch = 0; iPatch < nPatch; iPatch++) {
1748 
1749  patch = (AliEMCALTriggerPatchInfo*)fTriggerPatchInfo->At( iPatch );
1750  if (patch->IsMainTrigger()) {
1751  if(doSimpleOffline){
1752  if(patch->IsOfflineSimple()){
1753  switch(trigger){
1754  case kTriggerLevel0:
1755  // option not yet implemented in the trigger maker
1756  if(patch->IsLevel0()) selected = patch;
1757  break;
1758  case kTriggerLevel1Jet:
1759  if(patch->IsJetHighSimple() || patch->IsJetLowSimple()){
1760  if(!selected) selected = patch;
1761  else if(patch->GetADCOfflineAmp() > selected->GetADCOfflineAmp()) selected = patch;
1762  }
1763  break;
1764  case kTriggerLevel1Gamma:
1765  if(patch->IsGammaHighSimple() || patch->IsGammaLowSimple()){
1766  if(!selected) selected = patch;
1767  else if(patch->GetADCOfflineAmp() > selected->GetADCOfflineAmp()) selected = patch;
1768  }
1769  break;
1770  default: // Silence compiler warnings
1771  AliError("Untreated case: Main Patch is recalculated; should be in 'else' branch");
1772  };
1773  }
1774  } else { // Not OfflineSimple
1775  switch(trigger){
1776  case kTriggerLevel0:
1777  if(patch->IsLevel0()) selected = patch;
1778  break;
1779  case kTriggerLevel1Jet:
1780  if(patch->IsJetHigh() || patch->IsJetLow()){
1781  if(!selected) selected = patch;
1782  else if (patch->GetADCAmp() > selected->GetADCAmp())
1783  selected = patch;
1784  }
1785  break;
1786  case kTriggerLevel1Gamma:
1787  if(patch->IsGammaHigh() || patch->IsGammaLow()){
1788  if(!selected) selected = patch;
1789  else if (patch->GetADCAmp() > selected->GetADCAmp())
1790  selected = patch;
1791  }
1792  break;
1793  default:
1794  AliError("Untreated case: Main Patch is recalculated; should be in 'else' branch");
1795  };
1796  }
1797  }
1798  else if ((trigger == kTriggerRecalcJet && patch->IsRecalcJet()) ||
1799  (trigger == kTriggerRecalcGamma && patch->IsRecalcGamma())) { // recalculated patches
1800  if (doSimpleOffline && patch->IsOfflineSimple()) {
1801  if(!selected) selected = patch;
1802  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
1803  selected = patch;
1804  }
1805  else if (!doSimpleOffline && !patch->IsOfflineSimple()) {
1806  if(!selected) selected = patch;
1807  else if (patch->GetADCAmp() > selected->GetADCAmp())
1808  selected = patch;
1809  }
1810  }
1811  }
1812  return selected;
1813 }
1814 
1816 {
1817  if (!(InputEvent()->FindListObject(obj->GetName()))) {
1818  InputEvent()->AddObject(obj);
1819  }
1820  else {
1821  if (!attempt) {
1822  AliFatal(Form("%s: Container with name %s already present. Aborting", GetName(), obj->GetName()));
1823  }
1824  }
1825 }
1826 
1828 {
1829 
1830  if (!fGeom) {
1831  AliWarning(Form("%s - AliAnalysisTaskEmcal::IsTrackInEmcalAcceptance - Geometry is not available!", GetName()));
1832  return kFALSE;
1833  }
1834 
1835  Double_t minPhi = fGeom->GetArm1PhiMin() - edges;
1836  Double_t maxPhi = fGeom->GetArm1PhiMax() + edges;
1837 
1838  if (part->Phi() > minPhi && part->Phi() < maxPhi) {
1839  return kTRUE;
1840  }
1841  else {
1842  return kFALSE;
1843  }
1844 }
1845 
1847 {
1848  axis->SetBinLabel(1, "NullObject");
1849  axis->SetBinLabel(2, "Pt");
1850  axis->SetBinLabel(3, "Acceptance");
1851  axis->SetBinLabel(4, "MCLabel");
1852  axis->SetBinLabel(5, "BitMap");
1853  axis->SetBinLabel(6, "HF cut");
1854  axis->SetBinLabel(7, "Bit6");
1855  axis->SetBinLabel(8, "NotHybridTrack");
1856  axis->SetBinLabel(9, "MCFlag");
1857  axis->SetBinLabel(10, "MCGenerator");
1858  axis->SetBinLabel(11, "ChargeCut");
1859  axis->SetBinLabel(12, "MinDistanceTPCSectorEdge");
1860  axis->SetBinLabel(13, "Bit12");
1861  axis->SetBinLabel(14, "IsEMCal");
1862  axis->SetBinLabel(15, "Time");
1863  axis->SetBinLabel(16, "Energy");
1864  axis->SetBinLabel(17, "ExoticCut");
1865  axis->SetBinLabel(18, "Bit17");
1866  axis->SetBinLabel(19, "Area");
1867  axis->SetBinLabel(20, "AreaEmc");
1868  axis->SetBinLabel(21, "ZLeadingCh");
1869  axis->SetBinLabel(22, "ZLeadingEmc");
1870  axis->SetBinLabel(23, "NEF");
1871  axis->SetBinLabel(24, "MinLeadPt");
1872  axis->SetBinLabel(25, "MaxTrackPt");
1873  axis->SetBinLabel(26, "MaxClusterPt");
1874  axis->SetBinLabel(27, "Flavour");
1875  axis->SetBinLabel(28, "TagStatus");
1876  axis->SetBinLabel(29, "MinNConstituents");
1877  axis->SetBinLabel(30, "Bit29");
1878  axis->SetBinLabel(31, "Bit30");
1879  axis->SetBinLabel(32, "Bit31");
1880 }
1881 
1882 Double_t AliAnalysisTaskEmcal::GetParallelFraction(AliVParticle* part1, AliVParticle* part2)
1883 {
1884  TVector3 vect1(part1->Px(), part1->Py(), part1->Pz());
1885  TVector3 vect2(part2->Px(), part2->Py(), part2->Pz());
1886  Double_t z = (vect1 * vect2) / (vect2 * vect2);
1887  return z;
1888 }
1889 
1890 Double_t AliAnalysisTaskEmcal::GetParallelFraction(const TVector3& vect1, AliVParticle* part2)
1891 {
1892  TVector3 vect2(part2->Px(), part2->Py(), part2->Pz());
1893  Double_t z = (vect1 * vect2) / (vect2 * vect2);
1894  return z;
1895 }
1896 
1897 void AliAnalysisTaskEmcal::GetEtaPhiDiff(const AliVTrack *t, const AliVCluster *v, Double_t &phidiff, Double_t &etadiff)
1898 {
1899  phidiff = 999;
1900  etadiff = 999;
1901 
1902  if (!t||!v) return;
1903 
1904  Double_t veta = t->GetTrackEtaOnEMCal();
1905  Double_t vphi = t->GetTrackPhiOnEMCal();
1906 
1907  Float_t pos[3] = {0};
1908  v->GetPosition(pos);
1909  TVector3 cpos(pos);
1910  Double_t ceta = cpos.Eta();
1911  Double_t cphi = cpos.Phi();
1912  etadiff=veta-ceta;
1913  phidiff=TVector2::Phi_mpi_pi(vphi-cphi);
1914 }
1915 
1916 Byte_t AliAnalysisTaskEmcal::GetTrackType(const AliVTrack *t)
1917 {
1918  Byte_t ret = 0;
1919  if (t->TestBit(BIT(22)) && !t->TestBit(BIT(23)))
1920  ret = 1;
1921  else if (!t->TestBit(BIT(22)) && t->TestBit(BIT(23)))
1922  ret = 2;
1923  else if (t->TestBit(BIT(22)) && t->TestBit(BIT(23)))
1924  ret = 3;
1925  return ret;
1926 }
1927 
1928 Byte_t AliAnalysisTaskEmcal::GetTrackType(const AliAODTrack *aodTrack, UInt_t filterBit1, UInt_t filterBit2)
1929 {
1930 
1931  Int_t res = 0;
1932 
1933  if (aodTrack->TestFilterBit(filterBit1)) {
1934  res = 0;
1935  }
1936  else if (aodTrack->TestFilterBit(filterBit2)) {
1937  if ((aodTrack->GetStatus()&AliVTrack::kITSrefit)!=0) {
1938  res = 1;
1939  }
1940  else {
1941  res = 2;
1942  }
1943  }
1944  else {
1945  res = 3;
1946  }
1947 
1948  return res;
1949 }
1950 
1952 {
1953  if (!fPythiaInfo) {
1955  }
1956 
1957  AliStack* stack = mcEvent->Stack();
1958 
1959  const Int_t nprim = stack->GetNprimary();
1960  // reject if partons are missing from stack for some reason
1961  if (nprim < 8) return;
1962 
1963  TParticle *part6 = stack->Particle(6);
1964  TParticle *part7 = stack->Particle(7);
1965 
1966  fPythiaInfo->SetPartonFlag6(TMath::Abs(part6->GetPdgCode()));
1967  fPythiaInfo->SetParton6(part6->Pt(), part6->Eta(), part6->Phi(), part6->GetMass());
1968 
1969  fPythiaInfo->SetPartonFlag7(TMath::Abs(part7->GetPdgCode()));
1970  fPythiaInfo->SetParton7(part7->Pt(), part7->Eta(), part7->Phi(), part7->GetMass());
1971 
1972  AliGenPythiaEventHeader *pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(mcEvent->GenEventHeader());
1973  if(pythiaGenHeader){
1974  Float_t ptWeight=pythiaGenHeader->EventWeight();
1975  fPythiaInfo->SetPythiaEventWeight(ptWeight);}
1976 }
1977 
1979 {
1980  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1981  if (!mgr) {
1982  ::Error("AddAODHandler", "No analysis manager to connect to.");
1983  return NULL;
1984  }
1985 
1986  AliAODInputHandler* aodHandler = new AliAODInputHandler();
1987 
1988  AliVEventHandler *inputHandler=mgr->GetInputEventHandler();
1989  if (inputHandler && (inputHandler->IsA() == AliMultiInputEventHandler::Class())) {
1990  AliMultiInputEventHandler *multiInputHandler=(AliMultiInputEventHandler*)inputHandler;
1991  multiInputHandler->AddInputEventHandler(aodHandler);
1992  }
1993  else {
1994  if (!inputHandler) {
1995  mgr->SetInputEventHandler(aodHandler);
1996  }
1997  else {
1998  ::Error("AddAODHandler", "inputHandler is NOT null. AOD handler was NOT added !!!");
1999  return NULL;
2000  }
2001  }
2002 
2003  return aodHandler;
2004 }
2005 
2007 {
2008  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
2009  if (!mgr) {
2010  ::Error("AddESDHandler", "No analysis manager to connect to.");
2011  return NULL;
2012  }
2013 
2014  AliESDInputHandler *esdHandler = new AliESDInputHandler();
2015 
2016  AliVEventHandler *inputHandler=mgr->GetInputEventHandler();
2017  if (inputHandler && (inputHandler->IsA() == AliMultiInputEventHandler::Class())) {
2018  AliMultiInputEventHandler *multiInputHandler=(AliMultiInputEventHandler*)inputHandler;
2019  multiInputHandler->AddInputEventHandler(esdHandler);
2020  }
2021  else {
2022  if (!inputHandler) {
2023  mgr->SetInputEventHandler(esdHandler);
2024  }
2025  else {
2026  ::Error("AddESDHandler", "inputHandler is NOT null. ESD handler was NOT added !!!");
2027  return NULL;
2028  }
2029  }
2030 
2031  return esdHandler;
2032 }
2033 
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
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...
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