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