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