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