AliPhysics  eff0747 (eff0747)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
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 "AliAnalysisTaskEmcal.h"
17 
18 #include <TClonesArray.h>
19 #include <TList.h>
20 #include <TObject.h>
21 #include <TH1F.h>
22 #include <TProfile.h>
23 #include <TSystem.h>
24 #include <TFile.h>
25 #include <TChain.h>
26 #include <TKey.h>
27 
28 #include "AliStack.h"
29 #include "AliAODEvent.h"
30 #include "AliAnalysisManager.h"
31 #include "AliCentrality.h"
32 #include "AliEMCALGeometry.h"
33 #include "AliESDEvent.h"
34 #include "AliEmcalParticle.h"
35 #include "AliEventplane.h"
36 #include "AliInputEventHandler.h"
37 #include "AliLog.h"
38 #include "AliMCParticle.h"
39 #include "AliVCluster.h"
40 #include "AliVEventHandler.h"
41 #include "AliVParticle.h"
42 #include "AliAODTrack.h"
43 #include "AliVCaloTrigger.h"
44 #include "AliGenPythiaEventHeader.h"
45 #include "AliAODMCHeader.h"
46 #include "AliMCEvent.h"
47 #include "AliAnalysisUtils.h"
48 #include "AliEMCALTriggerPatchInfo.h"
49 #include "AliEmcalPythiaInfo.h"
50 
51 #include "AliMultSelection.h"
52 
54 
58 
63  AliAnalysisTaskSE("AliAnalysisTaskEmcal"),
64  fPythiaInfoName(""),
65  fForceBeamType(kNA),
66  fGeneralHistograms(kFALSE),
67  fInitialized(kFALSE),
68  fCreateHisto(kTRUE),
69  fCaloCellsName(),
70  fCaloTriggersName(),
71  fCaloTriggerPatchInfoName(),
72  fMinCent(-999),
73  fMaxCent(-999),
74  fMinVz(-999),
75  fMaxVz(999),
76  fTrackPtCut(0),
77  fMinNTrack(0),
78  fUseAliAnaUtils(kFALSE),
79  fRejectPileup(kFALSE),
80  fTklVsClusSPDCut(kFALSE),
81  fOffTrigger(AliVEvent::kAny),
82  fTrigClass(),
83  fTriggerTypeSel(kND),
84  fNbins(250),
85  fMinBinPt(0),
86  fMaxBinPt(250),
87  fMinPtTrackInEmcal(0),
88  fEventPlaneVsEmcal(-1),
89  fMinEventPlane(-1e6),
90  fMaxEventPlane(1e6),
91  fCentEst("V0M"),
92  fIsEmbedded(kFALSE),
93  fIsPythia(kFALSE),
94  fSelectPtHardBin(-999),
95  fMinMCLabel(0),
96  fMCLabelShift(0),
97  fNcentBins(4),
98  fNeedEmcalGeom(kTRUE),
99  fParticleCollArray(),
100  fClusterCollArray(),
101  fTriggers(0),
102  fEMCalTriggerMode(kOverlapWithLowThreshold),
103  fUseNewCentralityEstimation(kFALSE),
104  fGeneratePythiaInfoObject(kFALSE),
105  fAliAnalysisUtils(0x0),
106  fIsEsd(kFALSE),
107  fGeom(0),
108  fTracks(0),
109  fCaloClusters(0),
110  fCaloCells(0),
111  fCaloTriggers(0),
112  fTriggerPatchInfo(0),
113  fCent(0),
114  fCentBin(-1),
115  fEPV0(-1.0),
116  fEPV0A(-1.0),
117  fEPV0C(-1.0),
118  fNVertCont(0),
119  fBeamType(kNA),
120  fPythiaHeader(0),
121  fPtHard(0),
122  fPtHardBin(0),
123  fNTrials(0),
124  fXsection(0),
125  fPythiaInfo(0),
126  fOutput(0),
127  fHistEventCount(0),
128  fHistTrialsAfterSel(0),
129  fHistEventsAfterSel(0),
130  fHistXsectionAfterSel(0),
131  fHistTrials(0),
132  fHistEvents(0),
133  fHistXsection(0),
134  fHistPtHard(0),
135  fHistCentrality(0),
136  fHistZVertex(0),
137  fHistEventPlane(0),
138  fHistEventRejection(0),
139  fHistTriggerClasses(0)
140 {
141  fVertex[0] = 0;
142  fVertex[1] = 0;
143  fVertex[2] = 0;
144 
145  fParticleCollArray.SetOwner(kTRUE);
146  fClusterCollArray.SetOwner(kTRUE);
147 }
148 
154 AliAnalysisTaskEmcal::AliAnalysisTaskEmcal(const char *name, Bool_t histo) :
155  AliAnalysisTaskSE(name),
156  fPythiaInfoName(""),
157  fForceBeamType(kNA),
158  fGeneralHistograms(kFALSE),
159  fInitialized(kFALSE),
160  fCreateHisto(histo),
161  fCaloCellsName(),
162  fCaloTriggersName(),
163  fCaloTriggerPatchInfoName(),
164  fMinCent(-999),
165  fMaxCent(-999),
166  fMinVz(-999),
167  fMaxVz(999),
168  fTrackPtCut(0),
169  fMinNTrack(0),
170  fUseAliAnaUtils(kFALSE),
171  fRejectPileup(kFALSE),
172  fTklVsClusSPDCut(kFALSE),
173  fOffTrigger(AliVEvent::kAny),
174  fTrigClass(),
175  fTriggerTypeSel(kND),
176  fNbins(250),
177  fMinBinPt(0),
178  fMaxBinPt(250),
179  fMinPtTrackInEmcal(0),
180  fEventPlaneVsEmcal(-1),
181  fMinEventPlane(-1e6),
182  fMaxEventPlane(1e6),
183  fCentEst("V0M"),
184  fIsEmbedded(kFALSE),
185  fIsPythia(kFALSE),
186  fSelectPtHardBin(-999),
187  fMinMCLabel(0),
188  fMCLabelShift(0),
189  fNcentBins(4),
190  fNeedEmcalGeom(kTRUE),
191  fParticleCollArray(),
192  fClusterCollArray(),
193  fTriggers(0),
194  fEMCalTriggerMode(kOverlapWithLowThreshold),
195  fUseNewCentralityEstimation(kFALSE),
196  fGeneratePythiaInfoObject(kFALSE),
197  fAliAnalysisUtils(0x0),
198  fIsEsd(kFALSE),
199  fGeom(0),
200  fTracks(0),
201  fCaloClusters(0),
202  fCaloCells(0),
203  fCaloTriggers(0),
204  fTriggerPatchInfo(0),
205  fCent(0),
206  fCentBin(-1),
207  fEPV0(-1.0),
208  fEPV0A(-1.0),
209  fEPV0C(-1.0),
210  fNVertCont(0),
211  fBeamType(kNA),
212  fPythiaHeader(0),
213  fPtHard(0),
214  fPtHardBin(0),
215  fNTrials(0),
216  fXsection(0),
217  fPythiaInfo(0),
218  fOutput(0),
219  fHistEventCount(0),
220  fHistTrialsAfterSel(0),
221  fHistEventsAfterSel(0),
222  fHistXsectionAfterSel(0),
223  fHistTrials(0),
224  fHistEvents(0),
225  fHistXsection(0),
226  fHistPtHard(0),
227  fHistCentrality(0),
228  fHistZVertex(0),
229  fHistEventPlane(0),
230  fHistEventRejection(0),
231  fHistTriggerClasses(0)
232 {
233  fVertex[0] = 0;
234  fVertex[1] = 0;
235  fVertex[2] = 0;
236 
237  fParticleCollArray.SetOwner(kTRUE);
238  fClusterCollArray.SetOwner(kTRUE);
239 
240  if (fCreateHisto) {
241  DefineOutput(1, TList::Class());
242  }
243 }
244 
249 {
250 }
251 
252 void AliAnalysisTaskEmcal::SetClusPtCut(Double_t cut, Int_t c)
253 {
255  if (cont) cont->SetClusPtCut(cut);
256  else AliError(Form("%s in SetClusPtCut(...): container %d not found",GetName(),c));
257 }
258 
259 void AliAnalysisTaskEmcal::SetClusTimeCut(Double_t min, Double_t max, Int_t c)
260 {
262  if (cont) cont->SetClusTimeCut(min,max);
263  else AliError(Form("%s in SetClusTimeCut(...): container %d not found",GetName(),c));
264 }
265 
266 void AliAnalysisTaskEmcal::SetTrackPtCut(Double_t cut, Int_t c)
267 {
269  if (cont) cont->SetParticlePtCut(cut);
270  else AliError(Form("%s in SetTrackPtCut(...): container %d not found",GetName(),c));
271 
272  fTrackPtCut = cut;
273 }
274 
275 void AliAnalysisTaskEmcal::SetTrackEtaLimits(Double_t min, Double_t max, Int_t c)
276 {
278  if (cont) cont->SetParticleEtaLimits(min,max);
279  else AliError(Form("%s in SetTrackPtCut(...): container %d not found",GetName(),c));
280 }
281 
282 void AliAnalysisTaskEmcal::SetTrackPhiLimits(Double_t min, Double_t max, Int_t c)
283 {
285  if (cont) cont->SetParticlePhiLimits(min,max);
286  else AliError(Form("%s in SetTrackPhiLimits(...): container %d not found",GetName(),c));
287 }
288 
293 {
294  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
295  if (mgr) {
296  AliVEventHandler *evhand = mgr->GetInputEventHandler();
297  if (evhand) {
298  if (evhand->InheritsFrom("AliESDInputHandler")) {
299  fIsEsd = kTRUE;
300  }
301  else {
302  fIsEsd = kFALSE;
303  }
304  }
305  else {
306  AliError("Event handler not found!");
307  }
308  }
309  else {
310  AliError("Analysis manager not found!");
311  }
312 
313  if (!fCreateHisto)
314  return;
315 
316  OpenFile(1);
317  fOutput = new TList();
318  fOutput->SetOwner();
319 
320  if (fForceBeamType == kpp)
321  fNcentBins = 1;
322 
323  if (!fGeneralHistograms)
324  return;
325 
326  if (fIsPythia) {
327  fHistTrialsAfterSel = new TH1F("fHistTrialsAfterSel", "fHistTrialsAfterSel", 11, 0, 11);
328  fHistTrialsAfterSel->GetXaxis()->SetTitle("p_{T} hard bin");
329  fHistTrialsAfterSel->GetYaxis()->SetTitle("trials");
331 
332  fHistEventsAfterSel = new TH1F("fHistEventsAfterSel", "fHistEventsAfterSel", 11, 0, 11);
333  fHistEventsAfterSel->GetXaxis()->SetTitle("p_{T} hard bin");
334  fHistEventsAfterSel->GetYaxis()->SetTitle("total events");
336 
337  fHistXsectionAfterSel = new TProfile("fHistXsectionAfterSel", "fHistXsectionAfterSel", 11, 0, 11);
338  fHistXsectionAfterSel->GetXaxis()->SetTitle("p_{T} hard bin");
339  fHistXsectionAfterSel->GetYaxis()->SetTitle("xsection");
341 
342  fHistTrials = new TH1F("fHistTrials", "fHistTrials", 11, 0, 11);
343  fHistTrials->GetXaxis()->SetTitle("p_{T} hard bin");
344  fHistTrials->GetYaxis()->SetTitle("trials");
345  fOutput->Add(fHistTrials);
346 
347  fHistEvents = new TH1F("fHistEvents", "fHistEvents", 11, 0, 11);
348  fHistEvents->GetXaxis()->SetTitle("p_{T} hard bin");
349  fHistEvents->GetYaxis()->SetTitle("total events");
350  fOutput->Add(fHistEvents);
351 
352  fHistXsection = new TProfile("fHistXsection", "fHistXsection", 11, 0, 11);
353  fHistXsection->GetXaxis()->SetTitle("p_{T} hard bin");
354  fHistXsection->GetYaxis()->SetTitle("xsection");
355  fOutput->Add(fHistXsection);
356 
357  const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
358  const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
359 
360  for (Int_t i = 1; i < 12; i++) {
361  fHistTrialsAfterSel->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
362  fHistEventsAfterSel->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
363 
364  fHistTrials->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
365  fHistXsection->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
366  fHistEvents->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
367  }
368 
369  fHistPtHard = new TH1F("fHistPtHard", "fHistPtHard", fNbins*2, fMinBinPt, fMaxBinPt*4);
370  fHistPtHard->GetXaxis()->SetTitle("p_{T,hard} (GeV/c)");
371  fHistPtHard->GetYaxis()->SetTitle("counts");
372  fOutput->Add(fHistPtHard);
373  }
374 
375  fHistZVertex = new TH1F("fHistZVertex","Z vertex position", 60, -30, 30);
376  fHistZVertex->GetXaxis()->SetTitle("z");
377  fHistZVertex->GetYaxis()->SetTitle("counts");
378  fOutput->Add(fHistZVertex);
379 
380  if (fForceBeamType != kpp) {
381  fHistCentrality = new TH1F("fHistCentrality","Event centrality distribution", 200, 0, 100);
382  fHistCentrality->GetXaxis()->SetTitle("Centrality (%)");
383  fHistCentrality->GetYaxis()->SetTitle("counts");
384  fOutput->Add(fHistCentrality);
385 
386  fHistEventPlane = new TH1F("fHistEventPlane","Event plane", 120, -TMath::Pi(), TMath::Pi());
387  fHistEventPlane->GetXaxis()->SetTitle("event plane");
388  fHistEventPlane->GetYaxis()->SetTitle("counts");
389  fOutput->Add(fHistEventPlane);
390  }
391 
392  fHistEventRejection = new TH1F("fHistEventRejection","Reasons to reject event",20,0,20);
393 #if ROOT_VERSION_CODE < ROOT_VERSION(6,4,2)
394  fHistEventRejection->SetBit(TH1::kCanRebin);
395 #else
396  fHistEventRejection->SetCanExtend(TH1::kAllAxes);
397 #endif
398  fHistEventRejection->GetXaxis()->SetBinLabel(1,"PhysSel");
399  fHistEventRejection->GetXaxis()->SetBinLabel(2,"trigger");
400  fHistEventRejection->GetXaxis()->SetBinLabel(3,"trigTypeSel");
401  fHistEventRejection->GetXaxis()->SetBinLabel(4,"Cent");
402  fHistEventRejection->GetXaxis()->SetBinLabel(5,"vertex contr.");
403  fHistEventRejection->GetXaxis()->SetBinLabel(6,"Vz");
404  fHistEventRejection->GetXaxis()->SetBinLabel(7,"trackInEmcal");
405  fHistEventRejection->GetXaxis()->SetBinLabel(8,"minNTrack");
406  fHistEventRejection->GetXaxis()->SetBinLabel(9,"VtxSel2013pA");
407  fHistEventRejection->GetXaxis()->SetBinLabel(10,"PileUp");
408  fHistEventRejection->GetXaxis()->SetBinLabel(11,"EvtPlane");
409  fHistEventRejection->GetXaxis()->SetBinLabel(12,"SelPtHardBin");
410  fHistEventRejection->GetXaxis()->SetBinLabel(13,"Bkg evt");
411  fHistEventRejection->GetYaxis()->SetTitle("counts");
413 
414  fHistTriggerClasses = new TH1F("fHistTriggerClasses","fHistTriggerClasses",3,0,3);
415 #if ROOT_VERSION_CODE < ROOT_VERSION(6,4,2)
416  fHistTriggerClasses->SetBit(TH1::kCanRebin);
417 #else
418  fHistTriggerClasses->SetCanExtend(TH1::kAllAxes);
419 #endif
421 
422  fHistEventCount = new TH1F("fHistEventCount","fHistEventCount",2,0,2);
423  fHistEventCount->GetXaxis()->SetBinLabel(1,"Accepted");
424  fHistEventCount->GetXaxis()->SetBinLabel(2,"Rejected");
425  fHistEventCount->GetYaxis()->SetTitle("counts");
426  fOutput->Add(fHistEventCount);
427 
428  PostData(1, fOutput);
429 }
430 
432 {
433  if (fIsPythia) {
437  fHistPtHard->Fill(fPtHard);
438  }
439 
440  fHistZVertex->Fill(fVertex[2]);
441 
442  if (fForceBeamType != kpp) {
443  fHistCentrality->Fill(fCent);
444  fHistEventPlane->Fill(fEPV0);
445  }
446 
447  TObjArray* triggerClasses = InputEvent()->GetFiredTriggerClasses().Tokenize(" ");
448  TIter next(triggerClasses);
449  TObjString* triggerClass = 0;
450  while ((triggerClass = static_cast<TObjString*>(next()))) {
451  fHistTriggerClasses->Fill(triggerClass->GetString(), 1);
452  }
453  delete triggerClasses;
454  triggerClasses = 0;
455 
456  return kTRUE;
457 }
458 
464 {
465  if (!fInitialized)
466  ExecOnce();
467 
468  if (!fInitialized)
469  return;
470 
471  if (!RetrieveEventObjects())
472  return;
473 
474  if (IsEventSelected()) {
475  if (fGeneralHistograms) fHistEventCount->Fill("Accepted",1);
476  }
477  else {
478  if (fGeneralHistograms) fHistEventCount->Fill("Rejected",1);
479  return;
480  }
481 
483  if (!FillGeneralHistograms())
484  return;
485  }
486 
487  if (!Run())
488  return;
489 
490  if (fCreateHisto) {
491  if (!FillHistograms())
492  return;
493  }
494 
495  if (fCreateHisto && fOutput) {
496  // information for this iteration of the UserExec in the container
497  PostData(1, fOutput);
498  }
499 }
500 
509 Bool_t AliAnalysisTaskEmcal::AcceptCluster(AliVCluster *clus, Int_t c) const
510 {
511  AliWarning("AliAnalysisTaskEmcal::AcceptCluster method is deprecated. Please use GetCusterContainer(c)->AcceptCluster(clus).");
512 
513  if (!clus) return kFALSE;
514 
516  if (!cont) {
517  AliError(Form("%s:Container %d not found",GetName(),c));
518  return 0;
519  }
520 
521  return cont->AcceptCluster(clus);
522 }
523 
532 Bool_t AliAnalysisTaskEmcal::AcceptTrack(AliVParticle *track, Int_t c) const
533 {
534  AliWarning("AliAnalysisTaskEmcal::AcceptTrack method is deprecated. Please use GetParticleContainer(c)->AcceptParticle(clus).");
535 
536  if (!track) return kFALSE;
537 
539  if (!cont) {
540  AliError(Form("%s:Container %d not found",GetName(),c));
541  return 0;
542  }
543 
544  return cont->AcceptParticle(track);
545 }
546 
558 Bool_t AliAnalysisTaskEmcal::PythiaInfoFromFile(const char* currFile, Float_t &fXsec, Float_t &fTrials, Int_t &pthard)
559 {
560 
561  TString file(currFile);
562  fXsec = 0;
563  fTrials = 1;
564 
565  if (file.Contains(".zip#")) {
566  Ssiz_t pos1 = file.Index("root_archive",12,0,TString::kExact);
567  Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
568  Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
569  file.Replace(pos+1,pos2-pos1,"");
570  } else {
571  // not an archive take the basename....
572  file.ReplaceAll(gSystem->BaseName(file.Data()),"");
573  }
574  AliDebug(1,Form("File name: %s",file.Data()));
575 
576  // Get the pt hard bin
577  TString strPthard(file);
578 
579  strPthard.Remove(strPthard.Last('/'));
580  strPthard.Remove(strPthard.Last('/'));
581  if (strPthard.Contains("AOD")) strPthard.Remove(strPthard.Last('/'));
582  strPthard.Remove(0,strPthard.Last('/')+1);
583  if (strPthard.IsDec())
584  pthard = strPthard.Atoi();
585  else
586  AliWarning(Form("Could not extract file number from path %s", strPthard.Data()));
587 
588  // 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
589  TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root"));
590 
591  if (!fxsec) {
592  // next trial fetch the histgram file
593  fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
594  if (!fxsec) {
595  // not a severe condition but inciate that we have no information
596  return kFALSE;
597  } else {
598  // find the tlist we want to be independtent of the name so use the Tkey
599  TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
600  if (!key) {
601  fxsec->Close();
602  return kFALSE;
603  }
604  TList *list = dynamic_cast<TList*>(key->ReadObj());
605  if (!list) {
606  fxsec->Close();
607  return kFALSE;
608  }
609  fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
610  fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
611  fxsec->Close();
612  }
613  } else { // no tree pyxsec.root
614  TTree *xtree = (TTree*)fxsec->Get("Xsection");
615  if (!xtree) {
616  fxsec->Close();
617  return kFALSE;
618  }
619  UInt_t ntrials = 0;
620  Double_t xsection = 0;
621  xtree->SetBranchAddress("xsection",&xsection);
622  xtree->SetBranchAddress("ntrials",&ntrials);
623  xtree->GetEntry(0);
624  fTrials = ntrials;
625  fXsec = xsection;
626  fxsec->Close();
627  }
628  return kTRUE;
629 }
630 
636 {
638  return kTRUE;
639 
640  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
641  if (!tree) {
642  AliError(Form("%s - UserNotify: No current tree!",GetName()));
643  return kFALSE;
644  }
645 
646  Float_t xsection = 0;
647  Float_t trials = 0;
648  Int_t pthardbin = 0;
649 
650  TFile *curfile = tree->GetCurrentFile();
651  if (!curfile) {
652  AliError(Form("%s - UserNotify: No current file!",GetName()));
653  return kFALSE;
654  }
655 
656  TChain *chain = dynamic_cast<TChain*>(tree);
657  if (chain) tree = chain->GetTree();
658 
659  Int_t nevents = tree->GetEntriesFast();
660 
661  PythiaInfoFromFile(curfile->GetName(), xsection, trials, pthardbin);
662 
663  // TODO: Workaround
664  if ((pthardbin < 0) || (pthardbin > 10)) pthardbin = 0;
665 
666  fHistTrials->Fill(pthardbin, trials);
667  fHistXsection->Fill(pthardbin, xsection);
668  fHistEvents->Fill(pthardbin, nevents);
669 
670  return kTRUE;
671 }
672 
678 {
679  if (!fPythiaInfoName.IsNull() && !fPythiaInfo) {
680  fPythiaInfo = dynamic_cast<AliEmcalPythiaInfo*>(event->FindListObject(fPythiaInfoName));
681  if (!fPythiaInfo) {
682  AliError(Form("%s: Could not retrieve parton infos! %s!", GetName(), fPythiaInfoName.Data()));
683  return;
684  }
685  }
686 }
687 
692 {
693  if (!InputEvent()) {
694  AliError(Form("%s: Could not retrieve event! Returning!", GetName()));
695  return;
696  }
697 
698  LoadPythiaInfo(InputEvent());
699 
700  if (fNeedEmcalGeom) {
701  fGeom = AliEMCALGeometry::GetInstance();
702  if (!fGeom) {
703  AliError(Form("%s: Can not create geometry", GetName()));
704  return;
705  }
706  }
707 
708  if (fEventPlaneVsEmcal >= 0) {
709  if (fGeom) {
710  Double_t ep = (fGeom->GetArm1PhiMax() + fGeom->GetArm1PhiMin()) / 2 * TMath::DegToRad() + fEventPlaneVsEmcal - TMath::Pi();
711  fMinEventPlane = ep - TMath::Pi() / 4;
712  fMaxEventPlane = ep + TMath::Pi() / 4;
713  }
714  else {
715  AliWarning("Could not set event plane limits because EMCal geometry was not loaded!");
716  }
717  }
718 
719  //Load all requested track branches - each container knows name already
720  for (Int_t i =0; i<fParticleCollArray.GetEntriesFast(); i++) {
721  AliParticleContainer *cont = static_cast<AliParticleContainer*>(fParticleCollArray.At(i));
722  cont->SetArray(InputEvent());
723  }
724 
725  if (fParticleCollArray.GetEntriesFast()>0) {
727  if (!fTracks) {
728  AliError(Form("%s: Could not retrieve first track branch!", GetName()));
729  return;
730  }
731  }
732 
733  //Load all requested cluster branches - each container knows name already
734  for (Int_t i =0; i<fClusterCollArray.GetEntriesFast(); i++) {
735  AliClusterContainer *cont = static_cast<AliClusterContainer*>(fClusterCollArray.At(i));
736  cont->SetArray(InputEvent());
737  }
738 
739  if (fClusterCollArray.GetEntriesFast()>0) {
741  if (!fCaloClusters) {
742  AliError(Form("%s: Could not retrieve first cluster branch!", GetName()));
743  return;
744  }
745  }
746 
747  if (!fCaloCellsName.IsNull() && !fCaloCells) {
748  fCaloCells = dynamic_cast<AliVCaloCells*>(InputEvent()->FindListObject(fCaloCellsName));
749  if (!fCaloCells) {
750  AliError(Form("%s: Could not retrieve cells %s!", GetName(), fCaloCellsName.Data()));
751  return;
752  }
753  }
754 
755  if (!fCaloTriggersName.IsNull() && !fCaloTriggers) {
756  fCaloTriggers = dynamic_cast<AliVCaloTrigger*>(InputEvent()->FindListObject(fCaloTriggersName));
757  if (!fCaloTriggers) {
758  AliError(Form("%s: Could not retrieve calo triggers %s!", GetName(), fCaloTriggersName.Data()));
759  return;
760  }
761  }
762 
763  if (!fCaloTriggerPatchInfoName.IsNull() && !fTriggerPatchInfo) {
764  fTriggerPatchInfo = GetArrayFromEvent(fCaloTriggerPatchInfoName.Data(),"AliEMCALTriggerPatchInfo");
765  if (!fTriggerPatchInfo) {
766  AliError(Form("%s: Could not retrieve calo trigger patch info %s!", GetName(), fCaloTriggerPatchInfoName.Data()));
767  return;
768  }
769 
770  }
771 
772  fInitialized = kTRUE;
773 }
774 
781 {
782 
783  if (fForceBeamType != kNA)
784  return fForceBeamType;
785 
786  AliESDEvent *esd = dynamic_cast<AliESDEvent*>(InputEvent());
787  if (esd) {
788  const AliESDRun *run = esd->GetESDRun();
789  TString beamType = run->GetBeamType();
790  if (beamType == "p-p")
791  return kpp;
792  else if (beamType == "A-A")
793  return kAA;
794  else if (beamType == "p-A")
795  return kpA;
796  else
797  return kNA;
798  } else {
799  Int_t runNumber = InputEvent()->GetRunNumber();
800  if ((runNumber >= 136851 && runNumber <= 139517) || // LHC10h
801  (runNumber >= 166529 && runNumber <= 170593)) { // LHC11h
802  return kAA;
803  } else if ((runNumber>=188365 && runNumber <= 188366) || // LHC12g
804  (runNumber >= 195344 && runNumber <= 196608)) { // LHC13b-f
805  return kpA;
806  } else {
807  return kpp;
808  }
809  }
810 }
811 
819 {
820  if (!fTriggerPatchInfo)
821  return 0;
822 
823  //number of patches in event
824  Int_t nPatch = fTriggerPatchInfo->GetEntries();
825 
826  //loop over patches to define trigger type of event
827  Int_t nG1 = 0;
828  Int_t nG2 = 0;
829  Int_t nJ1 = 0;
830  Int_t nJ2 = 0;
831  Int_t nL0 = 0;
832  AliEMCALTriggerPatchInfo *patch;
833  for (Int_t iPatch = 0; iPatch < nPatch; iPatch++) {
834  patch = (AliEMCALTriggerPatchInfo*)fTriggerPatchInfo->At( iPatch );
835  if (patch->IsGammaHigh()) nG1++;
836  if (patch->IsGammaLow()) nG2++;
837  if (patch->IsJetHigh()) nJ1++;
838  if (patch->IsJetLow()) nJ2++;
839  if (patch->IsLevel0()) nL0++;
840  }
841 
842  AliDebug(2, "Patch summary: ");
843  AliDebug(2, Form("Number of patches: %d", nPatch));
844  AliDebug(2, Form("Jet: low[%d], high[%d]" ,nJ2, nJ1));
845  AliDebug(2, Form("Gamma: low[%d], high[%d]" ,nG2, nG1));
846 
847  ULong_t triggers(0);
848  if (nL0>0) SETBIT(triggers, kL0);
849  if (nG1>0) SETBIT(triggers, kG1);
850  if (nG2>0) SETBIT(triggers, kG2);
851  if (nJ1>0) SETBIT(triggers, kJ1);
852  if (nJ2>0) SETBIT(triggers, kJ2);
853  return triggers;
854 }
855 
862 {
863  //
864  if(trigger==kND) {
865  AliWarning(Form("%s: Requesting undefined trigger type!", GetName()));
866  return kFALSE;
867  }
868  //MV: removing this logic which as far as I can see doesn't make any sense
869  // if(trigger & kND){
870  // return fTriggers == 0;
871  // }
872  return TESTBIT(fTriggers, trigger);
873 }
874 
880 {
881  if (fOffTrigger != AliVEvent::kAny) {
882  UInt_t res = 0;
883  const AliESDEvent *eev = dynamic_cast<const AliESDEvent*>(InputEvent());
884  if (eev) {
885  res = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
886  } else {
887  const AliAODEvent *aev = dynamic_cast<const AliAODEvent*>(InputEvent());
888  if (aev) {
889  res = ((AliVAODHeader*)aev->GetHeader())->GetOfflineTrigger();
890  }
891  }
892  if ((res & fOffTrigger) == 0) {
893  if (fGeneralHistograms) fHistEventRejection->Fill("PhysSel",1);
894  return kFALSE;
895  }
896  }
897 
898  if (!fTrigClass.IsNull()) {
899  TString fired;
900  const AliESDEvent *eev = dynamic_cast<const AliESDEvent*>(InputEvent());
901  if (eev) {
902  fired = eev->GetFiredTriggerClasses();
903  } else {
904  const AliAODEvent *aev = dynamic_cast<const AliAODEvent*>(InputEvent());
905  if (aev) {
906  fired = aev->GetFiredTriggerClasses();
907  }
908  }
909  if (!fired.Contains("-B-")) {
910  if (fGeneralHistograms) fHistEventRejection->Fill("trigger",1);
911  return kFALSE;
912  }
913 
914  TObjArray *arr = fTrigClass.Tokenize("|");
915  if (!arr) {
916  if (fGeneralHistograms) fHistEventRejection->Fill("trigger",1);
917  return kFALSE;
918  }
919  Bool_t match = 0;
920  for (Int_t i=0;i<arr->GetEntriesFast();++i) {
921  TObject *obj = arr->At(i);
922  if (!obj)
923  continue;
924 
925  //Check if requested trigger was fired
926  TString objStr = obj->GetName();
928  (objStr.Contains("J1") || objStr.Contains("J2") || objStr.Contains("G1") || objStr.Contains("G2"))) {
929  // This is relevant for EMCal triggers with 2 thresholds
930  // If the kOverlapWithLowThreshold was requested than the overlap between the two triggers goes with the lower threshold trigger
931  TString trigType1 = "J1";
932  TString trigType2 = "J2";
933  if(objStr.Contains("G")) {
934  trigType1 = "G1";
935  trigType2 = "G2";
936  }
937  if(objStr.Contains(trigType2) && fired.Contains(trigType2.Data())) { //requesting low threshold + overlap
938  match = 1;
939  break;
940  }
941  else if(objStr.Contains(trigType1) && fired.Contains(trigType1.Data()) && !fired.Contains(trigType2.Data())) { //high threshold only
942  match = 1;
943  break;
944  }
945  }
946  else {
947  // If this is not an EMCal trigger, or no particular treatment of EMCal triggers was requested,
948  // simply check that the trigger was fired
949  if (fired.Contains(obj->GetName())) {
950  match = 1;
951  break;
952  }
953  }
954  }
955  delete arr;
956  if (!match) {
957  if (fGeneralHistograms) fHistEventRejection->Fill("trigger",1);
958  return kFALSE;
959  }
960  }
961 
962  if (fTriggerTypeSel != kND) {
964  if (fGeneralHistograms) fHistEventRejection->Fill("trigTypeSel",1);
965  return kFALSE;
966  }
967  }
968 
969  if ((fMinCent != -999) && (fMaxCent != -999)) {
970  if (fCent<fMinCent || fCent>fMaxCent) {
971  if (fGeneralHistograms) fHistEventRejection->Fill("Cent",1);
972  return kFALSE;
973  }
974  }
975 
976  if (fUseAliAnaUtils) {
977  if (!fAliAnalysisUtils)
978  fAliAnalysisUtils = new AliAnalysisUtils();
979  fAliAnalysisUtils->SetMinVtxContr(2);
980  fAliAnalysisUtils->SetMaxVtxZ(999);
981  if(fMinVz<-998.) fMinVz = -10.;
982  if(fMaxVz>998.) fMaxVz = 10.;
983 
984  if (!fAliAnalysisUtils->IsVertexSelected2013pA(InputEvent())) {
985  if (fGeneralHistograms) fHistEventRejection->Fill("VtxSel2013pA",1);
986  return kFALSE;
987  }
988 
989  if (fRejectPileup && fAliAnalysisUtils->IsPileUpEvent(InputEvent())) {
990  if (fGeneralHistograms) fHistEventRejection->Fill("PileUp",1);
991  return kFALSE;
992  }
993 
994  if(fTklVsClusSPDCut && fAliAnalysisUtils->IsSPDClusterVsTrackletBG(InputEvent())) {
995  if (fGeneralHistograms) fHistEventRejection->Fill("Bkg evt",1);
996  return kFALSE;
997  }
998  }
999 
1000  if ((fMinVz != -999) && (fMaxVz != -999)) {
1001  if (fNVertCont == 0 ) {
1002  if (fGeneralHistograms) fHistEventRejection->Fill("vertex contr.",1);
1003  return kFALSE;
1004  }
1005  Double_t vz = fVertex[2];
1006  if (vz<fMinVz || vz>fMaxVz) {
1007  if (fGeneralHistograms) fHistEventRejection->Fill("Vz",1);
1008  return kFALSE;
1009  }
1010  }
1011 
1012  if (fMinPtTrackInEmcal > 0 && fGeom) {
1013  Bool_t trackInEmcalOk = kFALSE;
1014  Int_t ntracks = GetNParticles(0);
1015  for (Int_t i = 0; i < ntracks; i++) {
1016  AliVParticle *track = GetAcceptParticleFromArray(i,0);
1017  if (!track)
1018  continue;
1019 
1020  Double_t phiMin = fGeom->GetArm1PhiMin() * TMath::DegToRad();
1021  Double_t phiMax = fGeom->GetArm1PhiMax() * TMath::DegToRad();
1022  Int_t runNumber = InputEvent()->GetRunNumber();
1023  if (runNumber>=177295 && runNumber<=197470) { //small SM masked in 2012 and 2013
1024  phiMin = 1.4;
1025  phiMax = TMath::Pi();
1026  }
1027 
1028  if (track->Eta() < fGeom->GetArm1EtaMin() || track->Eta() > fGeom->GetArm1EtaMax() || track->Phi() < phiMin || track->Phi() > phiMax)
1029  continue;
1030  if (track->Pt() > fMinPtTrackInEmcal) {
1031  trackInEmcalOk = kTRUE;
1032  break;
1033  }
1034  }
1035  if (!trackInEmcalOk) {
1036  if (fGeneralHistograms) fHistEventRejection->Fill("trackInEmcal",1);
1037  return kFALSE;
1038  }
1039  }
1040 
1041  if (fMinNTrack > 0) {
1042  Int_t nTracksAcc = 0;
1043  Int_t ntracks = GetNParticles(0);
1044  for (Int_t i = 0; i < ntracks; i++) {
1045  AliVParticle *track = GetAcceptParticleFromArray(i,0);
1046  if (!track)
1047  continue;
1048  if (track->Pt() > fTrackPtCut) {
1049  nTracksAcc++;
1050  if (nTracksAcc>=fMinNTrack)
1051  break;
1052  }
1053  }
1054  if (nTracksAcc<fMinNTrack) {
1055  if (fGeneralHistograms) fHistEventRejection->Fill("minNTrack",1);
1056  return kFALSE;
1057  }
1058  }
1059 
1060  if (!(fEPV0 > fMinEventPlane && fEPV0 <= fMaxEventPlane) &&
1061  !(fEPV0 + TMath::Pi() > fMinEventPlane && fEPV0 + TMath::Pi() <= fMaxEventPlane) &&
1062  !(fEPV0 - TMath::Pi() > fMinEventPlane && fEPV0 - TMath::Pi() <= fMaxEventPlane))
1063  {
1064  if (fGeneralHistograms) fHistEventRejection->Fill("EvtPlane",1);
1065  return kFALSE;
1066  }
1067 
1068  if (fSelectPtHardBin != -999 && fSelectPtHardBin != fPtHardBin) {
1069  if (fGeneralHistograms) fHistEventRejection->Fill("SelPtHardBin",1);
1070  return kFALSE;
1071  }
1072 
1073  return kTRUE;
1074 }
1075 
1082 TClonesArray *AliAnalysisTaskEmcal::GetArrayFromEvent(const char *name, const char *clname)
1083 {
1084  TClonesArray *arr = 0;
1085  TString sname(name);
1086  if (!sname.IsNull()) {
1087  arr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(sname));
1088  if (!arr) {
1089  AliWarning(Form("%s: Could not retrieve array with name %s!", GetName(), name));
1090  return 0;
1091  }
1092  } else {
1093  return 0;
1094  }
1095 
1096  if (!clname)
1097  return arr;
1098 
1099  TString objname(arr->GetClass()->GetName());
1100  TClass cls(objname);
1101  if (!cls.InheritsFrom(clname)) {
1102  AliWarning(Form("%s: Objects of type %s in %s are not inherited from %s!",
1103  GetName(), cls.GetName(), name, clname));
1104  return 0;
1105  }
1106  return arr;
1107 }
1108 
1114 {
1115  fVertex[0] = 0;
1116  fVertex[1] = 0;
1117  fVertex[2] = 0;
1118  fNVertCont = 0;
1119 
1120  if (fGeneratePythiaInfoObject && MCEvent()) {
1121  GeneratePythiaInfoObject(MCEvent());
1122  }
1123 
1124  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
1125  if (vert) {
1126  vert->GetXYZ(fVertex);
1127  fNVertCont = vert->GetNContributors();
1128  }
1129 
1130  fBeamType = GetBeamType();
1131 
1132  if (fBeamType == kAA || fBeamType == kpA ) {
1134  AliMultSelection *MultSelection = static_cast<AliMultSelection*>(InputEvent()->FindListObject("MultSelection"));
1135  if (MultSelection) {
1136  fCent = MultSelection->GetMultiplicityPercentile("V0M");
1137  }
1138  else {
1139  AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
1140  }
1141  }
1142  else { // old centrality estimation < 2015
1143  AliCentrality *aliCent = InputEvent()->GetCentrality();
1144  if (aliCent) {
1145  fCent = aliCent->GetCentralityPercentile(fCentEst.Data());
1146  }
1147  else {
1148  AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
1149  }
1150  }
1151 
1152  if (fNcentBins==4) {
1153  if (fCent >= 0 && fCent < 10) fCentBin = 0;
1154  else if (fCent >= 10 && fCent < 30) fCentBin = 1;
1155  else if (fCent >= 30 && fCent < 50) fCentBin = 2;
1156  else if (fCent >= 50 && fCent <= 100) fCentBin = 3;
1157  else {
1158  AliWarning(Form("%s: Negative centrality: %f. Assuming 99", GetName(), fCent));
1159  fCentBin = fNcentBins-1;
1160  }
1161  }
1162  else if (fNcentBins==5) { // for PbPb 2015
1163  if (fCent >= 0 && fCent < 10) fCentBin = 0;
1164  else if (fCent >= 10 && fCent < 30) fCentBin = 1;
1165  else if (fCent >= 30 && fCent < 50) fCentBin = 2;
1166  else if (fCent >= 50 && fCent <= 90) fCentBin = 3;
1167  else if (fCent > 90) {
1168  fCent = 99;
1169  fCentBin = 4;
1170  }
1171  else {
1172  AliWarning(Form("%s: Negative centrality: %f. Assuming 99", GetName(), fCent));
1173  fCentBin = fNcentBins-1;
1174  }
1175  }
1176  else {
1177  Double_t centWidth = (fMaxCent-fMinCent)/(Double_t)fNcentBins;
1178  if(centWidth>0.) {
1179  fCentBin = TMath::FloorNint(fCent/centWidth);
1180  }
1181  else {
1182  fCentBin = 0;
1183  }
1184  if (fCentBin>=fNcentBins) {
1185  AliWarning(Form("%s: fCentBin too large: cent = %f fCentBin = %d. Assuming 99", GetName(),fCent,fCentBin));
1186  fCentBin = fNcentBins-1;
1187  }
1188  }
1189 
1190  AliEventplane *aliEP = InputEvent()->GetEventplane();
1191  if (aliEP) {
1192  fEPV0 = aliEP->GetEventplane("V0" ,InputEvent());
1193  fEPV0A = aliEP->GetEventplane("V0A",InputEvent());
1194  fEPV0C = aliEP->GetEventplane("V0C",InputEvent());
1195  } else {
1196  AliWarning(Form("%s: Could not retrieve event plane information!", GetName()));
1197  }
1198  }
1199  else {
1200  fCent = 99;
1201  fCentBin = 0;
1202  }
1203 
1204  if (fIsPythia) {
1205 
1206  if (MCEvent()) {
1207  fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
1208  if (!fPythiaHeader) {
1209  // Check if AOD
1210  AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
1211 
1212  if (aodMCH) {
1213  for (UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++) {
1214  fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
1215  if (fPythiaHeader) break;
1216  }
1217  }
1218  }
1219  }
1220 
1221  if (fPythiaHeader) {
1222  fPtHard = fPythiaHeader->GetPtHard();
1223 
1224  const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
1225  const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
1226  for (fPtHardBin = 0; fPtHardBin < 11; fPtHardBin++) {
1227  if (fPtHard >= ptHardLo[fPtHardBin] && fPtHard < ptHardHi[fPtHardBin])
1228  break;
1229  }
1230 
1231  fXsection = fPythiaHeader->GetXsection();
1232  fNTrials = fPythiaHeader->Trials();
1233  }
1234  }
1235 
1237 
1238  AliEmcalContainer* cont = 0;
1239 
1240  TIter nextPartColl(&fParticleCollArray);
1241  while ((cont = static_cast<AliEmcalContainer*>(nextPartColl()))) cont->NextEvent();
1242 
1243  TIter nextClusColl(&fClusterCollArray);
1244  while ((cont = static_cast<AliParticleContainer*>(nextClusColl()))) cont->NextEvent();
1245 
1246  return kTRUE;
1247 }
1248 
1257 {
1258  if (TString(n).IsNull()) return 0;
1259 
1261 
1262  fParticleCollArray.Add(cont);
1263 
1264  return cont;
1265 }
1266 
1275 {
1276  if (TString(n).IsNull()) return 0;
1277 
1278  AliTrackContainer* cont = new AliTrackContainer(n);
1279 
1280  fParticleCollArray.Add(cont);
1281 
1282  return cont;
1283 }
1284 
1293 {
1294  if (TString(n).IsNull()) return 0;
1295 
1297 
1298  fParticleCollArray.Add(cont);
1299 
1300  return cont;
1301 }
1302 
1311 {
1312  if (TString(n).IsNull()) return 0;
1313 
1315 
1316  fClusterCollArray.Add(cont);
1317 
1318  return cont;
1319 }
1320 
1327 {
1328  if (i<0 || i>fParticleCollArray.GetEntriesFast()) return 0;
1329  AliParticleContainer *cont = static_cast<AliParticleContainer*>(fParticleCollArray.At(i));
1330  return cont;
1331 }
1332 
1339 {
1340  if (i<0 || i>fClusterCollArray.GetEntriesFast()) return 0;
1341  AliClusterContainer *cont = static_cast<AliClusterContainer*>(fClusterCollArray.At(i));
1342  return cont;
1343 }
1344 
1351 {
1352  AliParticleContainer *cont = static_cast<AliParticleContainer*>(fParticleCollArray.FindObject(name));
1353  return cont;
1354 }
1355 
1362 {
1363  AliClusterContainer *cont = static_cast<AliClusterContainer*>(fClusterCollArray.FindObject(name));
1364  return cont;
1365 }
1366 
1372 TClonesArray* AliAnalysisTaskEmcal::GetParticleArray(Int_t i) const
1373 {
1375  if (!cont) {
1376  AliError(Form("%s: Particle container %d not found",GetName(),i));
1377  return 0;
1378  }
1379  TString contName = cont->GetArrayName();
1380  return cont->GetArray();
1381 }
1382 
1388 TClonesArray* AliAnalysisTaskEmcal::GetClusterArray(Int_t i) const
1389 {
1391  if (!cont) {
1392  AliError(Form("%s:Cluster container %d not found",GetName(),i));
1393  return 0;
1394  }
1395  return cont->GetArray();
1396 }
1397 
1405 AliVParticle* AliAnalysisTaskEmcal::GetAcceptParticleFromArray(Int_t p, Int_t c) const
1406 {
1407 
1409  if (!cont) {
1410  AliError(Form("%s: Particle container %d not found",GetName(),c));
1411  return 0;
1412  }
1413  AliVParticle *vp = cont->GetAcceptParticle(p);
1414 
1415  return vp;
1416 }
1417 
1425 AliVCluster* AliAnalysisTaskEmcal::GetAcceptClusterFromArray(Int_t cl, Int_t c) const
1426 {
1428  if (!cont) {
1429  AliError(Form("%s: Cluster container %d not found",GetName(),c));
1430  return 0;
1431  }
1432  AliVCluster *vc = cont->GetAcceptCluster(cl);
1433 
1434  return vc;
1435 }
1436 
1443 {
1445  if (!cont) {
1446  AliError(Form("%s: Particle container %d not found",GetName(),i));
1447  return 0;
1448  }
1449  return cont->GetNEntries();
1450 }
1451 
1459 {
1461  if (!cont) {
1462  AliError(Form("%s: Cluster container %d not found",GetName(),i));
1463  return 0;
1464  }
1465  return cont->GetNEntries();
1466 }
1467 
1479 AliEMCALTriggerPatchInfo* AliAnalysisTaskEmcal::GetMainTriggerPatch(TriggerCategory trigger, Bool_t doSimpleOffline)
1480 {
1481 
1482  if (!fTriggerPatchInfo) {
1483  AliError(Form("%s: fTriggerPatchInfo not available",GetName()));
1484  return 0;
1485  }
1486 
1487  //number of patches in event
1488  Int_t nPatch = fTriggerPatchInfo->GetEntries();
1489 
1490  //extract main trigger patch(es)
1491  AliEMCALTriggerPatchInfo *patch(NULL), *selected(NULL);
1492  for (Int_t iPatch = 0; iPatch < nPatch; iPatch++) {
1493 
1494  patch = (AliEMCALTriggerPatchInfo*)fTriggerPatchInfo->At( iPatch );
1495  if (patch->IsMainTrigger()) {
1496  if(doSimpleOffline){
1497  if(patch->IsOfflineSimple()){
1498  switch(trigger){
1499  case kTriggerLevel0:
1500  // option not yet implemented in the trigger maker
1501  if(patch->IsLevel0()) selected = patch;
1502  break;
1503  case kTriggerLevel1Jet:
1504  if(patch->IsJetHighSimple() || patch->IsJetLowSimple()){
1505  if(!selected) selected = patch;
1506  else if(patch->GetADCOfflineAmp() > selected->GetADCOfflineAmp()) selected = patch;
1507  }
1508  break;
1509  case kTriggerLevel1Gamma:
1510  if(patch->IsGammaHighSimple() || patch->IsGammaLowSimple()){
1511  if(!selected) selected = patch;
1512  else if(patch->GetADCOfflineAmp() > selected->GetADCOfflineAmp()) selected = patch;
1513  }
1514  break;
1515  default: // Silence compiler warnings
1516  AliError("Untreated case: Main Patch is recalculated; should be in 'else' branch");
1517  };
1518  }
1519  } else { // Not OfflineSimple
1520  switch(trigger){
1521  case kTriggerLevel0:
1522  if(patch->IsLevel0()) selected = patch;
1523  break;
1524  case kTriggerLevel1Jet:
1525  if(patch->IsJetHigh() || patch->IsJetLow()){
1526  if(!selected) selected = patch;
1527  else if (patch->GetADCAmp() > selected->GetADCAmp())
1528  selected = patch;
1529  }
1530  break;
1531  case kTriggerLevel1Gamma:
1532  if(patch->IsGammaHigh() || patch->IsGammaLow()){
1533  if(!selected) selected = patch;
1534  else if (patch->GetADCAmp() > selected->GetADCAmp())
1535  selected = patch;
1536  }
1537  break;
1538  default:
1539  AliError("Untreated case: Main Patch is recalculated; should be in 'else' branch");
1540  };
1541  }
1542  }
1543  else if ((trigger == kTriggerRecalcJet && patch->IsRecalcJet()) ||
1544  (trigger == kTriggerRecalcGamma && patch->IsRecalcGamma())) { // recalculated patches
1545  if (doSimpleOffline && patch->IsOfflineSimple()) {
1546  if(!selected) selected = patch;
1547  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
1548  selected = patch;
1549  }
1550  else if (!doSimpleOffline && !patch->IsOfflineSimple()) {
1551  if(!selected) selected = patch;
1552  else if (patch->GetADCAmp() > selected->GetADCAmp())
1553  selected = patch;
1554  }
1555  }
1556  }
1557  return selected;
1558 }
1559 
1565 void AliAnalysisTaskEmcal::AddObjectToEvent(TObject *obj, Bool_t attempt)
1566 {
1567  if (!(InputEvent()->FindListObject(obj->GetName()))) {
1568  InputEvent()->AddObject(obj);
1569  }
1570  else {
1571  if (!attempt) {
1572  AliFatal(Form("%s: Container with name %s already present. Aborting", GetName(), obj->GetName()));
1573  }
1574  }
1575 }
1576 
1584 Bool_t AliAnalysisTaskEmcal::IsTrackInEmcalAcceptance(AliVParticle* part, Double_t edges) const
1585 {
1586 
1587  if (!fGeom) {
1588  AliWarning(Form("%s - AliAnalysisTaskEmcal::IsTrackInEmcalAcceptance - Geometry is not available!", GetName()));
1589  return kFALSE;
1590  }
1591 
1592  Double_t minPhi = fGeom->GetArm1PhiMin() - edges;
1593  Double_t maxPhi = fGeom->GetArm1PhiMax() + edges;
1594 
1595  if (part->Phi() > minPhi && part->Phi() < maxPhi) {
1596  return kTRUE;
1597  }
1598  else {
1599  return kFALSE;
1600  }
1601 }
1602 
1604 {
1605  axis->SetBinLabel(1, "NullObject");
1606  axis->SetBinLabel(2, "Pt");
1607  axis->SetBinLabel(3, "Acceptance");
1608  axis->SetBinLabel(4, "MCLabel");
1609  axis->SetBinLabel(5, "BitMap");
1610  axis->SetBinLabel(6, "HF cut");
1611  axis->SetBinLabel(7, "Bit6");
1612  axis->SetBinLabel(8, "NotHybridTrack");
1613  axis->SetBinLabel(9, "MCFlag");
1614  axis->SetBinLabel(10, "MCGenerator");
1615  axis->SetBinLabel(11, "ChargeCut");
1616  axis->SetBinLabel(12, "MinDistanceTPCSectorEdge");
1617  axis->SetBinLabel(13, "Bit12");
1618  axis->SetBinLabel(14, "IsEMCal");
1619  axis->SetBinLabel(15, "Time");
1620  axis->SetBinLabel(16, "Energy");
1621  axis->SetBinLabel(17, "ExoticCut");
1622  axis->SetBinLabel(18, "Bit17");
1623  axis->SetBinLabel(19, "Area");
1624  axis->SetBinLabel(20, "AreaEmc");
1625  axis->SetBinLabel(21, "ZLeadingCh");
1626  axis->SetBinLabel(22, "ZLeadingEmc");
1627  axis->SetBinLabel(23, "NEF");
1628  axis->SetBinLabel(24, "MinLeadPt");
1629  axis->SetBinLabel(25, "MaxTrackPt");
1630  axis->SetBinLabel(26, "MaxClusterPt");
1631  axis->SetBinLabel(27, "Flavour");
1632  axis->SetBinLabel(28, "TagStatus");
1633  axis->SetBinLabel(29, "MinNConstituents");
1634  axis->SetBinLabel(30, "Bit29");
1635  axis->SetBinLabel(31, "Bit30");
1636  axis->SetBinLabel(32, "Bit31");
1637 }
1638 
1645 Double_t AliAnalysisTaskEmcal::GetParallelFraction(AliVParticle* part1, AliVParticle* part2)
1646 {
1647  TVector3 vect1(part1->Px(), part1->Py(), part1->Pz());
1648  TVector3 vect2(part2->Px(), part2->Py(), part2->Pz());
1649  Double_t z = (vect1 * vect2) / (vect2 * vect2);
1650  return z;
1651 }
1652 
1659 Double_t AliAnalysisTaskEmcal::GetParallelFraction(const TVector3& vect1, AliVParticle* part2)
1660 {
1661  TVector3 vect2(part2->Px(), part2->Py(), part2->Pz());
1662  Double_t z = (vect1 * vect2) / (vect2 * vect2);
1663  return z;
1664 }
1665 
1674 void AliAnalysisTaskEmcal::GetEtaPhiDiff(const AliVTrack *t, const AliVCluster *v, Double_t &phidiff, Double_t &etadiff)
1675 {
1676  phidiff = 999;
1677  etadiff = 999;
1678 
1679  if (!t||!v) return;
1680 
1681  Double_t veta = t->GetTrackEtaOnEMCal();
1682  Double_t vphi = t->GetTrackPhiOnEMCal();
1683 
1684  Float_t pos[3] = {0};
1685  v->GetPosition(pos);
1686  TVector3 cpos(pos);
1687  Double_t ceta = cpos.Eta();
1688  Double_t cphi = cpos.Phi();
1689  etadiff=veta-ceta;
1690  phidiff=TVector2::Phi_mpi_pi(vphi-cphi);
1691 }
1692 
1698 Byte_t AliAnalysisTaskEmcal::GetTrackType(const AliVTrack *t)
1699 {
1700  Byte_t ret = 0;
1701  if (t->TestBit(BIT(22)) && !t->TestBit(BIT(23)))
1702  ret = 1;
1703  else if (!t->TestBit(BIT(22)) && t->TestBit(BIT(23)))
1704  ret = 2;
1705  else if (t->TestBit(BIT(22)) && t->TestBit(BIT(23)))
1706  ret = 3;
1707  return ret;
1708 }
1709 
1719 Byte_t AliAnalysisTaskEmcal::GetTrackType(const AliAODTrack *aodTrack, UInt_t filterBit1, UInt_t filterBit2)
1720 {
1721 
1722  Int_t res = 0;
1723 
1724  if (aodTrack->TestFilterBit(filterBit1)) {
1725  res = 0;
1726  }
1727  else if (aodTrack->TestFilterBit(filterBit2)) {
1728  if ((aodTrack->GetStatus()&AliVTrack::kITSrefit)!=0) {
1729  res = 1;
1730  }
1731  else {
1732  res = 2;
1733  }
1734  }
1735  else {
1736  res = 3;
1737  }
1738 
1739  return res;
1740 }
1741 
1747 {
1748  if (!fPythiaInfo) {
1750  }
1751 
1752  AliStack* stack = mcEvent->Stack();
1753 
1754  const Int_t nprim = stack->GetNprimary();
1755  // reject if partons are missing from stack for some reason
1756  if (nprim < 8) return;
1757 
1758  TParticle *part6 = stack->Particle(6);
1759  TParticle *part7 = stack->Particle(7);
1760 
1761  fPythiaInfo->SetPartonFlag6(TMath::Abs(part6->GetPdgCode()));
1762  fPythiaInfo->SetParton6(part6->Pt(), part6->Eta(), part6->Phi(), part6->GetMass());
1763 
1764  fPythiaInfo->SetPartonFlag7(TMath::Abs(part7->GetPdgCode()));
1765  fPythiaInfo->SetParton7(part7->Pt(), part7->Eta(), part7->Phi(), part7->GetMass());
1766 
1767  AliGenPythiaEventHeader *pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(mcEvent->GenEventHeader());
1768  if(pythiaGenHeader){
1769  Float_t ptWeight=pythiaGenHeader->EventWeight();
1770  fPythiaInfo->SetPythiaEventWeight(ptWeight);}
1771 
1772 
1773 }
ClassImp(AliAnalysisTaskTriggerRates) AliAnalysisTaskTriggerRates
Bool_t fGeneratePythiaInfoObject
Use new centrality estimation (for 2015 data)
TObjArray fClusterCollArray
particle/track collection array
void SetParticlePtCut(Double_t cut)
Bool_t fIsPythia
trigger, embedded signal
void SetParton7(Float_t pt, Float_t eta, Float_t phi, Float_t mass=0)
TH1 * fHistTrials
!trials from pyxsec.root
EMCAL Level1 gamma trigger, low threshold.
AliEmcalPythiaInfo * fPythiaInfo
!event parton info
Bool_t AcceptTrack(AliVParticle *track, Int_t c=0) const
EMCAL Level1 jet trigger, low threshold.
Bool_t HasTriggerType(TriggerType triggersel)
Int_t fNTrials
!event trials
UInt_t fOffTrigger
Apply tracklet-vs-cluster SPD cut to reject background events in pp.
Double_t fMinCent
trigger patch info array name
Double_t fTrackPtCut
max vertex for event selection
Recalculated jet trigger patch; does not need to be above trigger threshold.
Base task in the EMCAL framework.
void SetPartonFlag7(Int_t flag7)
TSystem * gSystem
Double_t fPtHard
!event pt hard
void SetTrackPtCut(Double_t cut, Int_t c=0)
static void GetEtaPhiDiff(const AliVTrack *t, const AliVCluster *v, Double_t &phidiff, Double_t &etadiff)
Double_t fMinBinPt
no. of pt bins
Double_t fEPV0
!event plane V0
TList * list
Bool_t fGeneralHistograms
forced beam type
virtual void SetArray(AliVEvent *event)
Bool_t AcceptCluster(AliVCluster *clus, Int_t c=0) const
Int_t fCentBin
!event centrality bin
TH1 * fHistEventsAfterSel
!total number of events per pt hard bin after selection
Double_t fMinPtTrackInEmcal
max pt in histograms
TH1 * fHistEventPlane
!event plane distribution
TList * fOutput
!output list
TH1 * fHistEvents
!total number of events per pt hard bin
void SetClusPtCut(Double_t cut, Int_t c=0)
AliClusterContainer * AddClusterContainer(const char *n)
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.
TObjArray fParticleCollArray
whether or not the task needs the emcal geometry
void SetTrackEtaLimits(Double_t min, Double_t max, Int_t c=0)
AliParticleContainer * GetParticleContainer(Int_t i=0) const
TProfile * fHistXsectionAfterSel
!x section from pythia header
TriggerType
Switch for EMCAL trigger types.
EMCalTriggerMode_t fEMCalTriggerMode
list of fired triggers
virtual Bool_t FillHistograms()
AliVCluster * GetAcceptCluster(Int_t i)
Int_t GetNParticles(Int_t i=0) const
TClonesArray * fCaloClusters
!clusters
Bool_t fUseNewCentralityEstimation
EMCal trigger selection mode.
Bool_t IsTrackInEmcalAcceptance(AliVParticle *part, Double_t edges=0.9) const
TH1 * fHistTriggerClasses
!number of events in each trigger class
Double_t fMaxVz
min vertex for event selection
void GeneratePythiaInfoObject(AliMCEvent *mcEvent)
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
name of calo triggers collection
TString fCaloTriggersName
name of calo cell collection
AliGenPythiaEventHeader * fPythiaHeader
!event Pythia header
void SetTrackPhiLimits(Double_t min, Double_t max, Int_t c=0)
AliParticleContainer * AddParticleContainer(const char *n)
AliAnalysisUtils * fAliAnalysisUtils
Generate Pythia info object.
BeamType fForceBeamType
name of pythia info object
Int_t fNcentBins
if MC label > fMCLabelShift, MC label -= fMCLabelShift
AliClusterContainer * GetClusterContainer(Int_t i=0) const
Base class for container structures within the EMCAL framework.
TriggerType fTriggerTypeSel
trigger class name for event selection
virtual Bool_t FillGeneralHistograms()
TString fTrigClass
offline trigger for event selection
TClonesArray * GetParticleArray(Int_t i=0) const
Double_t fMinVz
max centrality for event selection
BeamType fBeamType
!event beam type
Double_t fCent
!event centrality
TClonesArray * GetArray() const
Double_t fMinEventPlane
select events which have a certain event plane wrt the emcal
TString fCaloCellsName
whether or not create histograms
Int_t GetNClusters(Int_t i=0) const
Int_t fNVertCont
!event vertex number of contributors
EMCAL Level0 trigger.
EMCAL Level1 jet trigger, high threshold.
Int_t fSelectPtHardBin
trigger, if it is a PYTHIA production
AliMCParticleContainer * AddMCParticleContainer(const char *n)
static Double_t GetParallelFraction(AliVParticle *part1, AliVParticle *part2)
virtual Bool_t RetrieveEventObjects()
Bool_t fRejectPileup
used for LHC13* data: z-vtx, Ncontributors, z-vtx resolution cuts
TProfile * fHistXsection
!x section from pyxsec.root
Bool_t PythiaInfoFromFile(const char *currFile, Float_t &fXsec, Float_t &fTrials, Int_t &pthard)
void UserExec(Option_t *option)
void SetPartonFlag6(Int_t flag6)
AliVCaloCells * fCaloCells
!cells
const TString & GetArrayName() const
TClonesArray * GetArrayFromEvent(const char *name, const char *clname=0)
Double_t fEventPlaneVsEmcal
min pt track in emcal
virtual Bool_t IsEventSelected()
TH1 * fHistPtHard
!pt hard distribution
void SetParticleEtaLimits(Double_t min, Double_t max)
Double_t fMaxBinPt
min pt in histograms
Int_t fPtHardBin
!event pt hard bin
TClonesArray * fTracks
!tracks
TH1 * fHistTrialsAfterSel
!total number of trials per pt hard bin after selection
void LoadPythiaInfo(AliVEvent *event)
Bool_t fIsEsd
!whether it's an ESD analysis
Double_t fVertex[3]
!event vertex
AliTrackContainer * AddTrackContainer(const char *n)
Bool_t fCreateHisto
whether or not the task has been already initialized
Store some informaion about a Pythia eventThis class is used to store some information about a Pythia...
TFile * file
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
Double_t fEPV0A
!event plane V0A
virtual Bool_t AcceptCluster(Int_t i)
TString fCentEst
maximum event plane value
TString fPythiaInfoName
phi value used to distinguish between DCal and EMCal
Declaration of class AliEmcalPythiaInfo.
void SetClusPtCut(Double_t cut)
EMCAL Level1 gamma trigger, high threshold.
void AddObjectToEvent(TObject *obj, Bool_t attempt=kFALSE)
AliVCaloTrigger * fCaloTriggers
!calo triggers
void SetRejectionReasonLabels(TAxis *axis)
TH1 * fHistZVertex
!z vertex position
Int_t fMinNTrack
cut on track pt in event selection
Int_t GetNEntries() const
static Byte_t GetTrackType(const AliVTrack *t)
Bool_t fUseAliAnaUtils
minimum nr of tracks in event with pT>fTrackPtCut
void SetClusTimeCut(Double_t min, Double_t max, Int_t c=0)
ULong_t fTriggers
cluster collection array
Double_t fMaxEventPlane
minimum event plane value
void SetPythiaEventWeight(Float_t ptWeight)
Float_t fXsection
!x-section from pythia header
Bool_t fInitialized
whether or not it should fill some general histograms
TH1 * fHistEventCount
!incoming and selected events
virtual AliVParticle * GetAcceptParticle(Int_t i=-1)
Double_t fMaxCent
min 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
virtual Bool_t AcceptParticle(const AliVParticle *vp)
Container structure for EMCAL clusters.
virtual void NextEvent()
Container for MC-true particles within the EMCAL framework.
Bool_t fNeedEmcalGeom
how many centrality bins
AliVCluster * GetAcceptClusterFromArray(Int_t cl, Int_t c=0) const
Int_t fNbins
trigger type to select based on trigger patches
Bool_t fTklVsClusSPDCut
Reject pilup using function AliAnalysisUtils::IsPileUpEvent()
AliEMCALTriggerPatchInfo * GetMainTriggerPatch(TriggerCategory triggersel=kTriggerLevel1Jet, Bool_t doOfflinSimple=kFALSE)
static Double_t fgkEMCalDCalPhiDivide