AliPhysics  a0db429 (a0db429)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros
AliAnalysisTaskEmcal.cxx
Go to the documentation of this file.
1 //
2 // Emcal base analysis task.
3 //
4 // Author: S.Aiola, M. Verweij
5 
6 #include <RVersion.h>
7 #include "AliAnalysisTaskEmcal.h"
8 
9 #include <TClonesArray.h>
10 #include <TList.h>
11 #include <TObject.h>
12 #include <TH1F.h>
13 #include <TProfile.h>
14 #include <TSystem.h>
15 #include <TFile.h>
16 #include <TChain.h>
17 #include <TKey.h>
18 
19 #include "AliAODEvent.h"
20 #include "AliAnalysisManager.h"
21 #include "AliCentrality.h"
22 #include "AliEMCALGeometry.h"
23 #include "AliESDEvent.h"
24 #include "AliEmcalParticle.h"
25 #include "AliEventplane.h"
26 #include "AliInputEventHandler.h"
27 #include "AliLog.h"
28 #include "AliMCParticle.h"
29 #include "AliVCluster.h"
30 #include "AliVEventHandler.h"
31 #include "AliVParticle.h"
32 #include "AliAODTrack.h"
33 #include "AliVCaloTrigger.h"
34 #include "AliGenPythiaEventHeader.h"
35 #include "AliAODMCHeader.h"
36 #include "AliMCEvent.h"
37 #include "AliAnalysisUtils.h"
38 #include "AliEmcalTriggerPatchInfoAP.h"
39 
40 #include "AliParticleContainer.h"
41 #include "AliClusterContainer.h"
42 
43 #include "AliMultSelection.h"
44 
46 
48 
49 //________________________________________________________________________
51  AliAnalysisTaskSE("AliAnalysisTaskEmcal"),
52  fForceBeamType(kNA),
53  fGeneralHistograms(kFALSE),
54  fInitialized(kFALSE),
55  fCreateHisto(kTRUE),
56  fCaloCellsName(),
57  fCaloTriggersName(),
58  fCaloTriggerPatchInfoName(),
59  fMinCent(-999),
60  fMaxCent(-999),
61  fMinVz(-999),
62  fMaxVz(-999),
63  fTrackPtCut(0),
64  fMinNTrack(0),
65  fUseAliAnaUtils(kFALSE),
66  fRejectPileup(kFALSE),
67  fTklVsClusSPDCut(kFALSE),
68  fOffTrigger(AliVEvent::kAny),
69  fTrigClass(),
70  fTriggerTypeSel(kND),
71  fNbins(500),
72  fMinBinPt(0),
73  fMaxBinPt(250),
74  fMinPtTrackInEmcal(0),
75  fEventPlaneVsEmcal(-1),
76  fMinEventPlane(-1e6),
77  fMaxEventPlane(1e6),
78  fCentEst("V0M"),
79  fIsEmbedded(kFALSE),
80  fIsPythia(kFALSE),
81  fSelectPtHardBin(-999),
82  fMinMCLabel(0),
83  fMCLabelShift(0),
84  fNcentBins(4),
85  fNeedEmcalGeom(kTRUE),
86  fParticleCollArray(),
87  fClusterCollArray(),
88  fTriggers(0),
89  fEMCalTriggerMode(kOverlapWithLowThreshold),
90  fUseNewCentralityEstimation(kFALSE),
91  fAliAnalysisUtils(0x0),
92  fIsEsd(kFALSE),
93  fGeom(0),
94  fTracks(0),
95  fCaloClusters(0),
96  fCaloCells(0),
97  fCaloTriggers(0),
98  fTriggerPatchInfo(0),
99  fCent(0),
100  fCentBin(-1),
101  fEPV0(-1.0),
102  fEPV0A(-1.0),
103  fEPV0C(-1.0),
104  fNVertCont(0),
105  fBeamType(kNA),
106  fPythiaHeader(0),
107  fPtHard(0),
108  fPtHardBin(0),
109  fNTrials(0),
110  fXsection(0),
111  fOutput(0),
112  fHistEventCount(0),
113  fHistTrialsAfterSel(0),
114  fHistEventsAfterSel(0),
115  fHistXsectionAfterSel(0),
116  fHistTrials(0),
117  fHistEvents(0),
118  fHistXsection(0),
119  fHistPtHard(0),
120  fHistCentrality(0),
121  fHistZVertex(0),
122  fHistEventPlane(0),
123  fHistEventRejection(0),
124  fHistTriggerClasses(0)
125 {
126  // Default constructor.
127 
128  fVertex[0] = 0;
129  fVertex[1] = 0;
130  fVertex[2] = 0;
131 
132  fParticleCollArray.SetOwner(kTRUE);
133  fClusterCollArray.SetOwner(kTRUE);
134 }
135 
136 //________________________________________________________________________
137 AliAnalysisTaskEmcal::AliAnalysisTaskEmcal(const char *name, Bool_t histo) :
138  AliAnalysisTaskSE(name),
139  fForceBeamType(kNA),
140  fGeneralHistograms(kFALSE),
141  fInitialized(kFALSE),
142  fCreateHisto(histo),
143  fCaloCellsName(),
144  fCaloTriggersName(),
145  fCaloTriggerPatchInfoName(),
146  fMinCent(-999),
147  fMaxCent(-999),
148  fMinVz(-999),
149  fMaxVz(-999),
150  fTrackPtCut(0),
151  fMinNTrack(0),
152  fUseAliAnaUtils(kFALSE),
153  fRejectPileup(kFALSE),
154  fTklVsClusSPDCut(kFALSE),
155  fOffTrigger(AliVEvent::kAny),
156  fTrigClass(),
157  fTriggerTypeSel(kND),
158  fNbins(500),
159  fMinBinPt(0),
160  fMaxBinPt(250),
161  fMinPtTrackInEmcal(0),
162  fEventPlaneVsEmcal(-1),
163  fMinEventPlane(-1e6),
164  fMaxEventPlane(1e6),
165  fCentEst("V0M"),
166  fIsEmbedded(kFALSE),
167  fIsPythia(kFALSE),
168  fSelectPtHardBin(-999),
169  fMinMCLabel(0),
170  fMCLabelShift(0),
171  fNcentBins(4),
172  fNeedEmcalGeom(kTRUE),
173  fParticleCollArray(),
174  fClusterCollArray(),
175  fTriggers(0),
176  fEMCalTriggerMode(kOverlapWithLowThreshold),
177  fUseNewCentralityEstimation(kFALSE),
178  fAliAnalysisUtils(0x0),
179  fIsEsd(kFALSE),
180  fGeom(0),
181  fTracks(0),
182  fCaloClusters(0),
183  fCaloCells(0),
184  fCaloTriggers(0),
185  fTriggerPatchInfo(0),
186  fCent(0),
187  fCentBin(-1),
188  fEPV0(-1.0),
189  fEPV0A(-1.0),
190  fEPV0C(-1.0),
191  fNVertCont(0),
192  fBeamType(kNA),
193  fPythiaHeader(0),
194  fPtHard(0),
195  fPtHardBin(0),
196  fNTrials(0),
197  fXsection(0),
198  fOutput(0),
199  fHistEventCount(0),
200  fHistTrialsAfterSel(0),
201  fHistEventsAfterSel(0),
202  fHistXsectionAfterSel(0),
203  fHistTrials(0),
204  fHistEvents(0),
205  fHistXsection(0),
206  fHistPtHard(0),
207  fHistCentrality(0),
208  fHistZVertex(0),
209  fHistEventPlane(0),
210  fHistEventRejection(0),
211  fHistTriggerClasses(0)
212 {
213  // Standard constructor.
214 
215  fVertex[0] = 0;
216  fVertex[1] = 0;
217  fVertex[2] = 0;
218 
219  fParticleCollArray.SetOwner(kTRUE);
220  fClusterCollArray.SetOwner(kTRUE);
221 
222  if (fCreateHisto) {
223  DefineOutput(1, TList::Class());
224  }
225 }
226 
227 //________________________________________________________________________
229 {
230  // Destructor
231 }
232 
233 //________________________________________________________________________
234 void AliAnalysisTaskEmcal::SetClusPtCut(Double_t cut, Int_t c)
235 {
237  if (cont) cont->SetClusPtCut(cut);
238  else AliError(Form("%s in SetClusPtCut(...): container %d not found",GetName(),c));
239 }
240 
241 //________________________________________________________________________
242 void AliAnalysisTaskEmcal::SetClusTimeCut(Double_t min, Double_t max, Int_t c)
243 {
245  if (cont) cont->SetClusTimeCut(min,max);
246  else AliError(Form("%s in SetClusTimeCut(...): container %d not found",GetName(),c));
247 }
248 
249 //________________________________________________________________________
250 void AliAnalysisTaskEmcal::SetTrackPtCut(Double_t cut, Int_t c)
251 {
253  if (cont) cont->SetParticlePtCut(cut);
254  else AliError(Form("%s in SetTrackPtCut(...): container %d not found",GetName(),c));
255 
256  fTrackPtCut = cut;
257 }
258 
259 //________________________________________________________________________
260 void AliAnalysisTaskEmcal::SetTrackEtaLimits(Double_t min, Double_t max, Int_t c)
261 {
263  if (cont) cont->SetParticleEtaLimits(min,max);
264  else AliError(Form("%s in SetTrackPtCut(...): container %d not found",GetName(),c));
265 }
266 
267 //________________________________________________________________________
268 void AliAnalysisTaskEmcal::SetTrackPhiLimits(Double_t min, Double_t max, Int_t c)
269 {
271  if (cont) cont->SetParticlePhiLimits(min,max);
272  else AliError(Form("%s in SetTrackPhiLimits(...): container %d not found",GetName(),c));
273 }
274 
275 //________________________________________________________________________
277 {
278  // Create user output.
279 
280  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
281  if (mgr) {
282  AliVEventHandler *evhand = mgr->GetInputEventHandler();
283  if (evhand) {
284  if (evhand->InheritsFrom("AliESDInputHandler")) {
285  fIsEsd = kTRUE;
286  }
287  else {
288  fIsEsd = kFALSE;
289  }
290  }
291  else {
292  AliError("Event handler not found!");
293  }
294  }
295  else {
296  AliError("Analysis manager not found!");
297  }
298 
299  if (!fCreateHisto)
300  return;
301 
302  OpenFile(1);
303  fOutput = new TList();
304  fOutput->SetOwner();
305 
306  if (fForceBeamType == kpp)
307  fNcentBins = 1;
308 
309  if (!fGeneralHistograms)
310  return;
311 
312  if (fIsPythia) {
313  fHistTrialsAfterSel = new TH1F("fHistTrialsAfterSel", "fHistTrialsAfterSel", 11, 0, 11);
314  fHistTrialsAfterSel->GetXaxis()->SetTitle("p_{T} hard bin");
315  fHistTrialsAfterSel->GetYaxis()->SetTitle("trials");
317 
318  fHistEventsAfterSel = new TH1F("fHistEventsAfterSel", "fHistEventsAfterSel", 11, 0, 11);
319  fHistEventsAfterSel->GetXaxis()->SetTitle("p_{T} hard bin");
320  fHistEventsAfterSel->GetYaxis()->SetTitle("total events");
322 
323  fHistXsectionAfterSel = new TProfile("fHistXsectionAfterSel", "fHistXsectionAfterSel", 11, 0, 11);
324  fHistXsectionAfterSel->GetXaxis()->SetTitle("p_{T} hard bin");
325  fHistXsectionAfterSel->GetYaxis()->SetTitle("xsection");
327 
328  fHistTrials = new TH1F("fHistTrials", "fHistTrials", 11, 0, 11);
329  fHistTrials->GetXaxis()->SetTitle("p_{T} hard bin");
330  fHistTrials->GetYaxis()->SetTitle("trials");
331  fOutput->Add(fHistTrials);
332 
333  fHistEvents = new TH1F("fHistEvents", "fHistEvents", 11, 0, 11);
334  fHistEvents->GetXaxis()->SetTitle("p_{T} hard bin");
335  fHistEvents->GetYaxis()->SetTitle("total events");
336  fOutput->Add(fHistEvents);
337 
338  fHistXsection = new TProfile("fHistXsection", "fHistXsection", 11, 0, 11);
339  fHistXsection->GetXaxis()->SetTitle("p_{T} hard bin");
340  fHistXsection->GetYaxis()->SetTitle("xsection");
341  fOutput->Add(fHistXsection);
342 
343  const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
344  const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
345 
346  for (Int_t i = 1; i < 12; i++) {
347  fHistTrialsAfterSel->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
348  fHistEventsAfterSel->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
349 
350  fHistTrials->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
351  fHistXsection->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
352  fHistEvents->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
353  }
354 
355  fHistPtHard = new TH1F("fHistPtHard", "fHistPtHard", fNbins*2, fMinBinPt, fMaxBinPt*4);
356  fHistPtHard->GetXaxis()->SetTitle("p_{T,hard} (GeV/c)");
357  fHistPtHard->GetYaxis()->SetTitle("counts");
358  fOutput->Add(fHistPtHard);
359  }
360 
361  fHistZVertex = new TH1F("fHistZVertex","Z vertex position", 60, -30, 30);
362  fHistZVertex->GetXaxis()->SetTitle("z");
363  fHistZVertex->GetYaxis()->SetTitle("counts");
364  fOutput->Add(fHistZVertex);
365 
366  if (fForceBeamType != kpp) {
367  fHistCentrality = new TH1F("fHistCentrality","Event centrality distribution", 200, 0, 100);
368  fHistCentrality->GetXaxis()->SetTitle("Centrality (%)");
369  fHistCentrality->GetYaxis()->SetTitle("counts");
370  fOutput->Add(fHistCentrality);
371 
372  fHistEventPlane = new TH1F("fHistEventPlane","Event plane", 120, -TMath::Pi(), TMath::Pi());
373  fHistEventPlane->GetXaxis()->SetTitle("event plane");
374  fHistEventPlane->GetYaxis()->SetTitle("counts");
375  fOutput->Add(fHistEventPlane);
376  }
377 
378  fHistEventRejection = new TH1F("fHistEventRejection","Reasons to reject event",20,0,20);
379 #if ROOT_VERSION_CODE < ROOT_VERSION(6,4,2)
380  fHistEventRejection->SetBit(TH1::kCanRebin);
381 #else
382  fHistEventRejection->SetCanExtend(TH1::kAllAxes);
383 #endif
384  fHistEventRejection->GetXaxis()->SetBinLabel(1,"PhysSel");
385  fHistEventRejection->GetXaxis()->SetBinLabel(2,"trigger");
386  fHistEventRejection->GetXaxis()->SetBinLabel(3,"trigTypeSel");
387  fHistEventRejection->GetXaxis()->SetBinLabel(4,"Cent");
388  fHistEventRejection->GetXaxis()->SetBinLabel(5,"vertex contr.");
389  fHistEventRejection->GetXaxis()->SetBinLabel(6,"Vz");
390  fHistEventRejection->GetXaxis()->SetBinLabel(7,"trackInEmcal");
391  fHistEventRejection->GetXaxis()->SetBinLabel(8,"minNTrack");
392  fHistEventRejection->GetXaxis()->SetBinLabel(9,"VtxSel2013pA");
393  fHistEventRejection->GetXaxis()->SetBinLabel(10,"PileUp");
394  fHistEventRejection->GetXaxis()->SetBinLabel(11,"EvtPlane");
395  fHistEventRejection->GetXaxis()->SetBinLabel(12,"SelPtHardBin");
396  fHistEventRejection->GetXaxis()->SetBinLabel(13,"Bkg evt");
397  fHistEventRejection->GetYaxis()->SetTitle("counts");
399 
400  fHistTriggerClasses = new TH1F("fHistTriggerClasses","fHistTriggerClasses",3,0,3);
401 #if ROOT_VERSION_CODE < ROOT_VERSION(6,4,2)
402  fHistTriggerClasses->SetBit(TH1::kCanRebin);
403 #else
404  fHistTriggerClasses->SetCanExtend(TH1::kAllAxes);
405 #endif
407 
408  fHistEventCount = new TH1F("fHistEventCount","fHistEventCount",2,0,2);
409  fHistEventCount->GetXaxis()->SetBinLabel(1,"Accepted");
410  fHistEventCount->GetXaxis()->SetBinLabel(2,"Rejected");
411  fHistEventCount->GetYaxis()->SetTitle("counts");
412  fOutput->Add(fHistEventCount);
413 
414  PostData(1, fOutput);
415 }
416 
417 //________________________________________________________________________
419 {
420  if (fIsPythia) {
424  fHistPtHard->Fill(fPtHard);
425  }
426 
427  fHistZVertex->Fill(fVertex[2]);
428 
429  if (fForceBeamType != kpp) {
430  fHistCentrality->Fill(fCent);
431  fHistEventPlane->Fill(fEPV0);
432  }
433 
434  TObjArray* triggerClasses = InputEvent()->GetFiredTriggerClasses().Tokenize(" ");
435  TIter next(triggerClasses);
436  TObjString* triggerClass = 0;
437  while ((triggerClass = static_cast<TObjString*>(next()))) {
438  fHistTriggerClasses->Fill(triggerClass->GetString(), 1);
439  }
440  delete triggerClasses;
441  triggerClasses = 0;
442 
443  return kTRUE;
444 }
445 
446 //________________________________________________________________________
448 {
449  // Main loop, called for each event.
450 
451  if (!fInitialized)
452  ExecOnce();
453 
454  if (!fInitialized)
455  return;
456 
457  if (!RetrieveEventObjects())
458  return;
459 
460  if (IsEventSelected()) {
461  if (fGeneralHistograms) fHistEventCount->Fill("Accepted",1);
462  }
463  else {
464  if (fGeneralHistograms) fHistEventCount->Fill("Rejected",1);
465  return;
466  }
467 
469  if (!FillGeneralHistograms())
470  return;
471  }
472 
473  if (!Run())
474  return;
475 
476  if (fCreateHisto) {
477  if (!FillHistograms())
478  return;
479  }
480 
481  if (fCreateHisto && fOutput) {
482  // information for this iteration of the UserExec in the container
483  PostData(1, fOutput);
484  }
485 }
486 
487 //________________________________________________________________________
488 Bool_t AliAnalysisTaskEmcal::AcceptCluster(AliVCluster *clus, Int_t c) const
489 {
490  // Return true if cluster is accepted.
491 
492  if (!clus)
493  return kFALSE;
494 
496  if (!cont) {
497  AliError(Form("%s:Container %d not found",GetName(),c));
498  return 0;
499  }
500 
501  return cont->AcceptCluster(clus);
502 }
503 
504 //________________________________________________________________________
505 Bool_t AliAnalysisTaskEmcal::AcceptTrack(AliVParticle *track, Int_t c) const
506 {
507  // Return true if track is accepted.
508 
509  if (!track)
510  return kFALSE;
511 
513  if (!cont) {
514  AliError(Form("%s:Container %d not found",GetName(),c));
515  return 0;
516  }
517 
518  return cont->AcceptParticle(track);
519 }
520 
521 //________________________________________________________________________
522 Bool_t AliAnalysisTaskEmcal::PythiaInfoFromFile(const char* currFile, Float_t &fXsec, Float_t &fTrials, Int_t &pthard)
523 {
524  //
525  // Get the cross section and the trails either from pyxsec.root or from pysec_hists.root
526  // Get the pt hard bin from the file path
527  // This is to called in Notify and should provide the path to the AOD/ESD file
528  // (Partially copied from AliAnalysisHelperJetTasks)
529 
530  TString file(currFile);
531  fXsec = 0;
532  fTrials = 1;
533 
534  if (file.Contains(".zip#")) {
535  Ssiz_t pos1 = file.Index("root_archive",12,0,TString::kExact);
536  Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
537  Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
538  file.Replace(pos+1,pos2-pos1,"");
539  } else {
540  // not an archive take the basename....
541  file.ReplaceAll(gSystem->BaseName(file.Data()),"");
542  }
543  AliDebug(1,Form("File name: %s",file.Data()));
544 
545  // Get the pt hard bin
546  TString strPthard(file);
547 
548  strPthard.Remove(strPthard.Last('/'));
549  strPthard.Remove(strPthard.Last('/'));
550  if (strPthard.Contains("AOD")) strPthard.Remove(strPthard.Last('/'));
551  strPthard.Remove(0,strPthard.Last('/')+1);
552  if (strPthard.IsDec())
553  pthard = strPthard.Atoi();
554  else
555  AliWarning(Form("Could not extract file number from path %s", strPthard.Data()));
556 
557  // 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
558  TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root"));
559 
560  if (!fxsec) {
561  // next trial fetch the histgram file
562  fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
563  if (!fxsec) {
564  // not a severe condition but inciate that we have no information
565  return kFALSE;
566  } else {
567  // find the tlist we want to be independtent of the name so use the Tkey
568  TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
569  if (!key) {
570  fxsec->Close();
571  return kFALSE;
572  }
573  TList *list = dynamic_cast<TList*>(key->ReadObj());
574  if (!list) {
575  fxsec->Close();
576  return kFALSE;
577  }
578  fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
579  fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
580  fxsec->Close();
581  }
582  } else { // no tree pyxsec.root
583  TTree *xtree = (TTree*)fxsec->Get("Xsection");
584  if (!xtree) {
585  fxsec->Close();
586  return kFALSE;
587  }
588  UInt_t ntrials = 0;
589  Double_t xsection = 0;
590  xtree->SetBranchAddress("xsection",&xsection);
591  xtree->SetBranchAddress("ntrials",&ntrials);
592  xtree->GetEntry(0);
593  fTrials = ntrials;
594  fXsec = xsection;
595  fxsec->Close();
596  }
597  return kTRUE;
598 }
599 
600 //________________________________________________________________________
602 {
603  // Called when file changes.
604 
606  return kTRUE;
607 
608  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
609  if (!tree) {
610  AliError(Form("%s - UserNotify: No current tree!",GetName()));
611  return kFALSE;
612  }
613 
614  Float_t xsection = 0;
615  Float_t trials = 0;
616  Int_t pthardbin = 0;
617 
618  TFile *curfile = tree->GetCurrentFile();
619  if (!curfile) {
620  AliError(Form("%s - UserNotify: No current file!",GetName()));
621  return kFALSE;
622  }
623 
624  TChain *chain = dynamic_cast<TChain*>(tree);
625  if (chain) tree = chain->GetTree();
626 
627  Int_t nevents = tree->GetEntriesFast();
628 
629  PythiaInfoFromFile(curfile->GetName(), xsection, trials, pthardbin);
630 
631  // TODO: Workaround
632  if ((pthardbin < 0) || (pthardbin > 10)) pthardbin = 0;
633 
634  fHistTrials->Fill(pthardbin, trials);
635  fHistXsection->Fill(pthardbin, xsection);
636  fHistEvents->Fill(pthardbin, nevents);
637 
638  return kTRUE;
639 }
640 
641 //________________________________________________________________________
643 {
644  // Init the analysis.
645 
646  if (!InputEvent()) {
647  AliError(Form("%s: Could not retrieve event! Returning!", GetName()));
648  return;
649  }
650 
651  if (fNeedEmcalGeom) {
652  fGeom = AliEMCALGeometry::GetInstance();
653  if (!fGeom) {
654  AliError(Form("%s: Can not create geometry", GetName()));
655  return;
656  }
657  }
658 
659  if (fEventPlaneVsEmcal >= 0) {
660  if (fGeom) {
661  Double_t ep = (fGeom->GetArm1PhiMax() + fGeom->GetArm1PhiMin()) / 2 * TMath::DegToRad() + fEventPlaneVsEmcal - TMath::Pi();
662  fMinEventPlane = ep - TMath::Pi() / 4;
663  fMaxEventPlane = ep + TMath::Pi() / 4;
664  }
665  else {
666  AliWarning("Could not set event plane limits because EMCal geometry was not loaded!");
667  }
668  }
669 
670  //Load all requested track branches - each container knows name already
671  for (Int_t i =0; i<fParticleCollArray.GetEntriesFast(); i++) {
672  AliParticleContainer *cont = static_cast<AliParticleContainer*>(fParticleCollArray.At(i));
673  cont->SetArray(InputEvent());
674  }
675 
676  if (fParticleCollArray.GetEntriesFast()>0) {
678  if (!fTracks) {
679  AliError(Form("%s: Could not retrieve first track branch!", GetName()));
680  return;
681  }
682  }
683 
684  //Load all requested cluster branches - each container knows name already
685  for (Int_t i =0; i<fClusterCollArray.GetEntriesFast(); i++) {
686  AliClusterContainer *cont = static_cast<AliClusterContainer*>(fClusterCollArray.At(i));
687  cont->SetArray(InputEvent());
688  }
689 
690  if (fClusterCollArray.GetEntriesFast()>0) {
692  if (!fCaloClusters) {
693  AliError(Form("%s: Could not retrieve first cluster branch!", GetName()));
694  return;
695  }
696  }
697 
698  if (!fCaloCellsName.IsNull() && !fCaloCells) {
699  fCaloCells = dynamic_cast<AliVCaloCells*>(InputEvent()->FindListObject(fCaloCellsName));
700  if (!fCaloCells) {
701  AliError(Form("%s: Could not retrieve cells %s!", GetName(), fCaloCellsName.Data()));
702  return;
703  }
704  }
705 
706  if (!fCaloTriggersName.IsNull() && !fCaloTriggers) {
707  fCaloTriggers = dynamic_cast<AliVCaloTrigger*>(InputEvent()->FindListObject(fCaloTriggersName));
708  if (!fCaloTriggers) {
709  AliError(Form("%s: Could not retrieve calo triggers %s!", GetName(), fCaloTriggersName.Data()));
710  return;
711  }
712  }
713 
714  if (!fCaloTriggerPatchInfoName.IsNull() && !fTriggerPatchInfo) {
715  fTriggerPatchInfo = GetArrayFromEvent(fCaloTriggerPatchInfoName.Data(),"AliEmcalTriggerPatchInfo");
716  if (!fTriggerPatchInfo) {
717  AliError(Form("%s: Could not retrieve calo trigger patch info %s!", GetName(), fCaloTriggerPatchInfoName.Data()));
718  return;
719  }
720 
721  }
722 
723  fInitialized = kTRUE;
724 }
725 
726 //_____________________________________________________
728 {
729  // Get beam type : pp-AA-pA
730  // ESDs have it directly, AODs get it from hardcoded run number ranges
731 
732  if (fForceBeamType != kNA)
733  return fForceBeamType;
734 
735  AliESDEvent *esd = dynamic_cast<AliESDEvent*>(InputEvent());
736  if (esd) {
737  const AliESDRun *run = esd->GetESDRun();
738  TString beamType = run->GetBeamType();
739  if (beamType == "p-p")
740  return kpp;
741  else if (beamType == "A-A")
742  return kAA;
743  else if (beamType == "p-A")
744  return kpA;
745  else
746  return kNA;
747  } else {
748  Int_t runNumber = InputEvent()->GetRunNumber();
749  if ((runNumber >= 136851 && runNumber <= 139517) || // LHC10h
750  (runNumber >= 166529 && runNumber <= 170593)) { // LHC11h
751  return kAA;
752  } else if ((runNumber>=188365 && runNumber <= 188366) || // LHC12g
753  (runNumber >= 195344 && runNumber <= 196608)) { // LHC13b-f
754  return kpA;
755  } else {
756  return kpp;
757  }
758  }
759 }
760 
761 //________________________________________________________________________
763 {
764  if (!fTriggerPatchInfo)
765  return 0;
766 
767  //number of patches in event
768  Int_t nPatch = fTriggerPatchInfo->GetEntries();
769 
770  //loop over patches to define trigger type of event
771  Int_t nG1 = 0;
772  Int_t nG2 = 0;
773  Int_t nJ1 = 0;
774  Int_t nJ2 = 0;
775  Int_t nL0 = 0;
776  AliEmcalTriggerPatchInfo *patch;
777  for (Int_t iPatch = 0; iPatch < nPatch; iPatch++) {
778  patch = (AliEmcalTriggerPatchInfo*)fTriggerPatchInfo->At( iPatch );
779  if (patch->IsGammaHigh()) nG1++;
780  if (patch->IsGammaLow()) nG2++;
781  if (patch->IsJetHigh()) nJ1++;
782  if (patch->IsJetLow()) nJ2++;
783  if (patch->IsLevel0()) nL0++;
784  }
785 
786  AliDebug(2, "Patch summary: ");
787  AliDebug(2, Form("Number of patches: %d", nPatch));
788  AliDebug(2, Form("Jet: low[%d], high[%d]" ,nJ2, nJ1));
789  AliDebug(2, Form("Gamma: low[%d], high[%d]" ,nG2, nG1));
790 
791  ULong_t triggers(0);
792  if (nL0>0) SETBIT(triggers, kL0);
793  if (nG1>0) SETBIT(triggers, kG1);
794  if (nG2>0) SETBIT(triggers, kG2);
795  if (nJ1>0) SETBIT(triggers, kJ1);
796  if (nJ2>0) SETBIT(triggers, kJ2);
797  return triggers;
798 }
799 
800 //________________________________________________________________________
802 {
803  // Check if event has a given trigger type
804  if(trigger==kND) {
805  AliWarning(Form("%s: Requesting undefined trigger type!", GetName()));
806  return kFALSE;
807  }
808  //MV: removing this logic which as far as I can see doesn't make any sense
809  // if(trigger & kND){
810  // return fTriggers == 0;
811  // }
812  return TESTBIT(fTriggers, trigger);
813 }
814 
815 //________________________________________________________________________
817 {
818  // Check if event is selected
819 
820  if (fOffTrigger != AliVEvent::kAny) {
821  UInt_t res = 0;
822  const AliESDEvent *eev = dynamic_cast<const AliESDEvent*>(InputEvent());
823  if (eev) {
824  res = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
825  } else {
826  const AliAODEvent *aev = dynamic_cast<const AliAODEvent*>(InputEvent());
827  if (aev) {
828  res = ((AliVAODHeader*)aev->GetHeader())->GetOfflineTrigger();
829  }
830  }
831  if ((res & fOffTrigger) == 0) {
832  if (fGeneralHistograms) fHistEventRejection->Fill("PhysSel",1);
833  return kFALSE;
834  }
835  }
836 
837  if (!fTrigClass.IsNull()) {
838  TString fired;
839  const AliESDEvent *eev = dynamic_cast<const AliESDEvent*>(InputEvent());
840  if (eev) {
841  fired = eev->GetFiredTriggerClasses();
842  } else {
843  const AliAODEvent *aev = dynamic_cast<const AliAODEvent*>(InputEvent());
844  if (aev) {
845  fired = aev->GetFiredTriggerClasses();
846  }
847  }
848  if (!fired.Contains("-B-")) {
849  if (fGeneralHistograms) fHistEventRejection->Fill("trigger",1);
850  return kFALSE;
851  }
852 
853  TObjArray *arr = fTrigClass.Tokenize("|");
854  if (!arr) {
855  if (fGeneralHistograms) fHistEventRejection->Fill("trigger",1);
856  return kFALSE;
857  }
858  Bool_t match = 0;
859  for (Int_t i=0;i<arr->GetEntriesFast();++i) {
860  TObject *obj = arr->At(i);
861  if (!obj)
862  continue;
863 
864  //Check if requested trigger was fired
865  TString objStr = obj->GetName();
867  (objStr.Contains("J1") || objStr.Contains("J2") || objStr.Contains("G1") || objStr.Contains("G2"))) {
868  // This is relevant for EMCal triggers with 2 thresholds
869  // If the kOverlapWithLowThreshold was requested than the overlap between the two triggers goes with the lower threshold trigger
870  TString trigType1 = "J1";
871  TString trigType2 = "J2";
872  if(objStr.Contains("G")) {
873  trigType1 = "G1";
874  trigType2 = "G2";
875  }
876  if(objStr.Contains(trigType2) && fired.Contains(trigType2.Data())) { //requesting low threshold + overlap
877  match = 1;
878  break;
879  }
880  else if(objStr.Contains(trigType1) && fired.Contains(trigType1.Data()) && !fired.Contains(trigType2.Data())) { //high threshold only
881  match = 1;
882  break;
883  }
884  }
885  else {
886  // If this is not an EMCal trigger, or no particular treatment of EMCal triggers was requested,
887  // simply check that the trigger was fired
888  if (fired.Contains(obj->GetName())) {
889  match = 1;
890  break;
891  }
892  }
893  }
894  delete arr;
895  if (!match) {
896  if (fGeneralHistograms) fHistEventRejection->Fill("trigger",1);
897  return kFALSE;
898  }
899  }
900 
901  if (fTriggerTypeSel != kND) {
903  if (fGeneralHistograms) fHistEventRejection->Fill("trigTypeSel",1);
904  return kFALSE;
905  }
906  }
907 
908  if ((fMinCent != -999) && (fMaxCent != -999)) {
909  if (fCent<fMinCent || fCent>fMaxCent) {
910  if (fGeneralHistograms) fHistEventRejection->Fill("Cent",1);
911  return kFALSE;
912  }
913  }
914 
915  if (fUseAliAnaUtils) {
916  if (!fAliAnalysisUtils)
917  fAliAnalysisUtils = new AliAnalysisUtils();
918  fAliAnalysisUtils->SetMinVtxContr(2);
919  fAliAnalysisUtils->SetMaxVtxZ(999);
920  if(fMinVz<-10.) fMinVz = -10.;
921  if(fMinVz>10.) fMaxVz = 10.;
922 
923  if (!fAliAnalysisUtils->IsVertexSelected2013pA(InputEvent())) {
924  if (fGeneralHistograms) fHistEventRejection->Fill("VtxSel2013pA",1);
925  return kFALSE;
926  }
927 
928  if (fRejectPileup && fAliAnalysisUtils->IsPileUpEvent(InputEvent())) {
929  if (fGeneralHistograms) fHistEventRejection->Fill("PileUp",1);
930  return kFALSE;
931  }
932 
933  if(fTklVsClusSPDCut && fAliAnalysisUtils->IsSPDClusterVsTrackletBG(InputEvent())) {
934  if (fGeneralHistograms) fHistEventRejection->Fill("Bkg evt",1);
935  return kFALSE;
936  }
937  }
938 
939  if ((fMinVz != -999) && (fMaxVz != -999)) {
940  if (fNVertCont == 0 ) {
941  if (fGeneralHistograms) fHistEventRejection->Fill("vertex contr.",1);
942  return kFALSE;
943  }
944  Double_t vz = fVertex[2];
945  if (vz<fMinVz || vz>fMaxVz) {
946  if (fGeneralHistograms) fHistEventRejection->Fill("Vz",1);
947  return kFALSE;
948  }
949  }
950 
951  if (fMinPtTrackInEmcal > 0 && fGeom) {
952  Bool_t trackInEmcalOk = kFALSE;
953  Int_t ntracks = GetNParticles(0);
954  for (Int_t i = 0; i < ntracks; i++) {
955  AliVParticle *track = GetAcceptParticleFromArray(i,0);
956  if (!track)
957  continue;
958 
959  Double_t phiMin = fGeom->GetArm1PhiMin() * TMath::DegToRad();
960  Double_t phiMax = fGeom->GetArm1PhiMax() * TMath::DegToRad();
961  Int_t runNumber = InputEvent()->GetRunNumber();
962  if (runNumber>=177295 && runNumber<=197470) { //small SM masked in 2012 and 2013
963  phiMin = 1.4;
964  phiMax = TMath::Pi();
965  }
966 
967  if (track->Eta() < fGeom->GetArm1EtaMin() || track->Eta() > fGeom->GetArm1EtaMax() || track->Phi() < phiMin || track->Phi() > phiMax)
968  continue;
969  if (track->Pt() > fMinPtTrackInEmcal) {
970  trackInEmcalOk = kTRUE;
971  break;
972  }
973  }
974  if (!trackInEmcalOk) {
975  if (fGeneralHistograms) fHistEventRejection->Fill("trackInEmcal",1);
976  return kFALSE;
977  }
978  }
979 
980  if (fMinNTrack > 0) {
981  Int_t nTracksAcc = 0;
982  Int_t ntracks = GetNParticles(0);
983  for (Int_t i = 0; i < ntracks; i++) {
984  AliVParticle *track = GetAcceptParticleFromArray(i,0);
985  if (!track)
986  continue;
987  if (track->Pt() > fTrackPtCut) {
988  nTracksAcc++;
989  if (nTracksAcc>=fMinNTrack)
990  break;
991  }
992  }
993  if (nTracksAcc<fMinNTrack) {
994  if (fGeneralHistograms) fHistEventRejection->Fill("minNTrack",1);
995  return kFALSE;
996  }
997  }
998 
999  if (!(fEPV0 > fMinEventPlane && fEPV0 <= fMaxEventPlane) &&
1000  !(fEPV0 + TMath::Pi() > fMinEventPlane && fEPV0 + TMath::Pi() <= fMaxEventPlane) &&
1001  !(fEPV0 - TMath::Pi() > fMinEventPlane && fEPV0 - TMath::Pi() <= fMaxEventPlane))
1002  {
1003  if (fGeneralHistograms) fHistEventRejection->Fill("EvtPlane",1);
1004  return kFALSE;
1005  }
1006 
1007  if (fSelectPtHardBin != -999 && fSelectPtHardBin != fPtHardBin) {
1008  if (fGeneralHistograms) fHistEventRejection->Fill("SelPtHardBin",1);
1009  return kFALSE;
1010  }
1011 
1012  return kTRUE;
1013 }
1014 
1015 //________________________________________________________________________
1016 TClonesArray *AliAnalysisTaskEmcal::GetArrayFromEvent(const char *name, const char *clname)
1017 {
1018  // Get array from event.
1019 
1020  TClonesArray *arr = 0;
1021  TString sname(name);
1022  if (!sname.IsNull()) {
1023  arr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(sname));
1024  if (!arr) {
1025  AliWarning(Form("%s: Could not retrieve array with name %s!", GetName(), name));
1026  return 0;
1027  }
1028  } else {
1029  return 0;
1030  }
1031 
1032  if (!clname)
1033  return arr;
1034 
1035  TString objname(arr->GetClass()->GetName());
1036  TClass cls(objname);
1037  if (!cls.InheritsFrom(clname)) {
1038  AliWarning(Form("%s: Objects of type %s in %s are not inherited from %s!",
1039  GetName(), cls.GetName(), name, clname));
1040  return 0;
1041  }
1042  return arr;
1043 }
1044 
1045 //________________________________________________________________________
1047 {
1048  // Retrieve objects from event.
1049 
1050  fVertex[0] = 0;
1051  fVertex[1] = 0;
1052  fVertex[2] = 0;
1053  fNVertCont = 0;
1054 
1055  const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
1056  if (vert) {
1057  vert->GetXYZ(fVertex);
1058  fNVertCont = vert->GetNContributors();
1059  }
1060 
1061  fBeamType = GetBeamType();
1062 
1063  if (fBeamType == kAA || fBeamType == kpA ) {
1065  AliMultSelection *MultSelection = static_cast<AliMultSelection*>(InputEvent()->FindListObject("MultSelection"));
1066  if (MultSelection) {
1067  fCent = MultSelection->GetMultiplicityPercentile("V0M");
1068  }
1069  else {
1070  AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
1071  }
1072  }
1073  else { // old centrality estimation < 2015
1074  AliCentrality *aliCent = InputEvent()->GetCentrality();
1075  if (aliCent) {
1076  fCent = aliCent->GetCentralityPercentile(fCentEst.Data());
1077  }
1078  else {
1079  AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
1080  }
1081  }
1082 
1083  if (fNcentBins==4) {
1084  if (fCent >= 0 && fCent < 10) fCentBin = 0;
1085  else if (fCent >= 10 && fCent < 30) fCentBin = 1;
1086  else if (fCent >= 30 && fCent < 50) fCentBin = 2;
1087  else if (fCent >= 50 && fCent <= 100) fCentBin = 3;
1088  else {
1089  AliWarning(Form("%s: Negative centrality: %f. Assuming 99", GetName(), fCent));
1090  fCentBin = fNcentBins-1;
1091  }
1092  }
1093  else if (fNcentBins==5) { // for PbPb 2015
1094  if (fCent >= 0 && fCent < 10) fCentBin = 0;
1095  else if (fCent >= 10 && fCent < 30) fCentBin = 1;
1096  else if (fCent >= 30 && fCent < 50) fCentBin = 2;
1097  else if (fCent >= 50 && fCent <= 90) fCentBin = 3;
1098  else if (fCent > 90) {
1099  fCent = 99;
1100  fCentBin = 4;
1101  }
1102  else {
1103  AliWarning(Form("%s: Negative centrality: %f. Assuming 99", GetName(), fCent));
1104  fCentBin = fNcentBins-1;
1105  }
1106  }
1107  else {
1108  Double_t centWidth = (fMaxCent-fMinCent)/(Double_t)fNcentBins;
1109  if(centWidth>0.) {
1110  fCentBin = TMath::FloorNint(fCent/centWidth);
1111  }
1112  else {
1113  fCentBin = 0;
1114  }
1115  if (fCentBin>=fNcentBins) {
1116  AliWarning(Form("%s: fCentBin too large: cent = %f fCentBin = %d. Assuming 99", GetName(),fCent,fCentBin));
1117  fCentBin = fNcentBins-1;
1118  }
1119  }
1120 
1121  AliEventplane *aliEP = InputEvent()->GetEventplane();
1122  if (aliEP) {
1123  fEPV0 = aliEP->GetEventplane("V0" ,InputEvent());
1124  fEPV0A = aliEP->GetEventplane("V0A",InputEvent());
1125  fEPV0C = aliEP->GetEventplane("V0C",InputEvent());
1126  } else {
1127  AliWarning(Form("%s: Could not retrieve event plane information!", GetName()));
1128  }
1129  }
1130 
1131  if (fIsPythia) {
1132 
1133  if (MCEvent()) {
1134  fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
1135  if (!fPythiaHeader) {
1136  // Check if AOD
1137  AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
1138 
1139  if (aodMCH) {
1140  for (UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++) {
1141  fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
1142  if (fPythiaHeader) break;
1143  }
1144  }
1145  }
1146  }
1147 
1148  if (fPythiaHeader) {
1149  fPtHard = fPythiaHeader->GetPtHard();
1150 
1151  const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
1152  const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
1153  for (fPtHardBin = 0; fPtHardBin < 11; fPtHardBin++) {
1154  if (fPtHard >= ptHardLo[fPtHardBin] && fPtHard < ptHardHi[fPtHardBin])
1155  break;
1156  }
1157 
1158  fXsection = fPythiaHeader->GetXsection();
1159  fNTrials = fPythiaHeader->Trials();
1160  }
1161  }
1162 
1164 
1165  return kTRUE;
1166 }
1167 
1168 //________________________________________________________________________
1170 {
1171  // Add particle container
1172  // will be called in AddTask macro
1173 
1174  TString tmp = TString(n);
1175  if (tmp.IsNull()) return 0;
1176 
1177  AliParticleContainer *cont = 0x0;
1178  cont = new AliParticleContainer();
1179  cont->SetArrayName(n);
1180  TString contName = cont->GetArrayName();
1181 
1182  fParticleCollArray.Add(cont);
1183 
1184  return cont;
1185 }
1186 
1187 //________________________________________________________________________
1189 {
1190  // Add cluster container
1191  // will be called in AddTask macro
1192 
1193  TString tmp = TString(n);
1194  if (tmp.IsNull()) return 0;
1195 
1196  AliClusterContainer *cont = 0x0;
1197  cont = new AliClusterContainer();
1198  cont->SetArrayName(n);
1199 
1200  fClusterCollArray.Add(cont);
1201 
1202  return cont;
1203 }
1204 
1205 //________________________________________________________________________
1207 {
1208  // Get i^th particle container
1209 
1210  if (i<0 || i>fParticleCollArray.GetEntriesFast()) return 0;
1211  AliParticleContainer *cont = static_cast<AliParticleContainer*>(fParticleCollArray.At(i));
1212  return cont;
1213 }
1214 
1215 //________________________________________________________________________
1217 {
1218  // Get i^th cluster container
1219 
1220  if (i<0 || i>fClusterCollArray.GetEntriesFast()) return 0;
1221  AliClusterContainer *cont = static_cast<AliClusterContainer*>(fClusterCollArray.At(i));
1222  return cont;
1223 }
1224 
1225 //________________________________________________________________________
1227 {
1228  // Get particle container with name
1229 
1230  AliParticleContainer *cont = static_cast<AliParticleContainer*>(fParticleCollArray.FindObject(name));
1231  return cont;
1232 }
1233 
1234 //________________________________________________________________________
1236 {
1237  // Get cluster container with name
1238 
1239  AliClusterContainer *cont = static_cast<AliClusterContainer*>(fClusterCollArray.FindObject(name));
1240  return cont;
1241 }
1242 
1243 //________________________________________________________________________
1244 TClonesArray* AliAnalysisTaskEmcal::GetParticleArray(Int_t i) const
1245 {
1246  // Get i^th TClonesArray with AliVParticle
1247 
1249  if (!cont) {
1250  AliError(Form("%s: Particle container %d not found",GetName(),i));
1251  return 0;
1252  }
1253  TString contName = cont->GetArrayName();
1254  return cont->GetArray();
1255 }
1256 
1257 //________________________________________________________________________
1258 TClonesArray* AliAnalysisTaskEmcal::GetClusterArray(Int_t i) const
1259 {
1260  // Get i^th TClonesArray with AliVCluster
1261 
1263  if (!cont) {
1264  AliError(Form("%s:Cluster container %d not found",GetName(),i));
1265  return 0;
1266  }
1267  return cont->GetArray();
1268 }
1269 
1270 //________________________________________________________________________
1271 AliVParticle* AliAnalysisTaskEmcal::GetAcceptParticleFromArray(Int_t p, Int_t c) const
1272 {
1273  // Get particle p if accepted from container c
1274  // If particle not accepted return 0
1275 
1277  if (!cont) {
1278  AliError(Form("%s: Particle container %d not found",GetName(),c));
1279  return 0;
1280  }
1281  AliVParticle *vp = cont->GetAcceptParticle(p);
1282 
1283  return vp;
1284 }
1285 
1286 //________________________________________________________________________
1287 AliVCluster* AliAnalysisTaskEmcal::GetAcceptClusterFromArray(Int_t cl, Int_t c) const
1288 {
1289  // Get particle p if accepted from container c
1290  // If particle not accepted return 0
1291 
1293  if (!cont) {
1294  AliError(Form("%s: Cluster container %d not found",GetName(),c));
1295  return 0;
1296  }
1297  AliVCluster *vc = cont->GetAcceptCluster(cl);
1298 
1299  return vc;
1300 }
1301 
1302 //________________________________________________________________________
1304 {
1305  // Get number of entries in particle array i
1306 
1308  if (!cont) {
1309  AliError(Form("%s: Particle container %d not found",GetName(),i));
1310  return 0;
1311  }
1312  return cont->GetNEntries();
1313 }
1314 
1315 //________________________________________________________________________
1317 {
1318  // Get number of entries in cluster array i
1319 
1321  if (!cont) {
1322  AliError(Form("%s: Cluster container %d not found",GetName(),i));
1323  return 0;
1324  }
1325  return cont->GetNEntries();
1326 }
1327 
1328 //________________________________________________________________________
1329 AliEmcalTriggerPatchInfo* AliAnalysisTaskEmcal::GetMainTriggerPatch(TriggerCategory trigger, Bool_t doSimpleOffline)
1330 {
1331  // Get main trigger match
1332  //
1333  // For the selection of the main trigger patch, high and low threshold triggers of a given category are grouped
1334  // If there are more than 1 main patch of a given trigger category (i.e. different high and low threshold patches),
1335  // the highest one according to the ADC value is taken. In case doSimpleOffline is true, then only the patches from
1336  // the simple offline trigger are used.
1337 
1338  if (!fTriggerPatchInfo) {
1339  AliError(Form("%s: fTriggerPatchInfo not available",GetName()));
1340  return 0;
1341  }
1342 
1343  //number of patches in event
1344  Int_t nPatch = fTriggerPatchInfo->GetEntries();
1345 
1346  //extract main trigger patch(es)
1347  AliEmcalTriggerPatchInfo *patch(NULL), *selected(NULL);
1348  for (Int_t iPatch = 0; iPatch < nPatch; iPatch++) {
1349 
1350  patch = (AliEmcalTriggerPatchInfo*)fTriggerPatchInfo->At( iPatch );
1351  if (patch->IsMainTrigger()) {
1352  if(doSimpleOffline){
1353  if(patch->IsOfflineSimple()){
1354  switch(trigger){
1355  case kTriggerLevel0:
1356  // option not yet implemented in the trigger maker
1357  if(patch->IsLevel0()) selected = patch;
1358  break;
1359  case kTriggerLevel1Jet:
1360  if(patch->IsJetHighSimple() || patch->IsJetLowSimple()){
1361  if(!selected) selected = patch;
1362  else if(patch->GetADCOfflineAmp() > selected->GetADCOfflineAmp()) selected = patch;
1363  }
1364  break;
1365  case kTriggerLevel1Gamma:
1366  if(patch->IsGammaHighSimple() || patch->IsGammaLowSimple()){
1367  if(!selected) selected = patch;
1368  else if(patch->GetADCOfflineAmp() > selected->GetADCOfflineAmp()) selected = patch;
1369  }
1370  break;
1371  default: // Silence compiler warnings
1372  AliError("Untreated case: Main Patch is recalculated; should be in 'else' branch");
1373  };
1374  }
1375  } else { // Not OfflineSimple
1376  switch(trigger){
1377  case kTriggerLevel0:
1378  if(patch->IsLevel0()) selected = patch;
1379  break;
1380  case kTriggerLevel1Jet:
1381  if(patch->IsJetHigh() || patch->IsJetLow()){
1382  if(!selected) selected = patch;
1383  else if (patch->GetADCAmp() > selected->GetADCAmp())
1384  selected = patch;
1385  }
1386  break;
1387  case kTriggerLevel1Gamma:
1388  if(patch->IsGammaHigh() || patch->IsGammaLow()){
1389  if(!selected) selected = patch;
1390  else if (patch->GetADCAmp() > selected->GetADCAmp())
1391  selected = patch;
1392  }
1393  break;
1394  default:
1395  AliError("Untreated case: Main Patch is recalculated; should be in 'else' branch");
1396  };
1397  }
1398  }
1399  else if ((trigger == kTriggerRecalcJet && patch->IsRecalcJet()) ||
1400  (trigger == kTriggerRecalcGamma && patch->IsRecalcGamma())) { // recalculated patches
1401  if (doSimpleOffline && patch->IsOfflineSimple()) {
1402  if(!selected) selected = patch;
1403  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
1404  selected = patch;
1405  }
1406  else if (!doSimpleOffline && !patch->IsOfflineSimple()) {
1407  if(!selected) selected = patch;
1408  else if (patch->GetADCAmp() > selected->GetADCAmp())
1409  selected = patch;
1410  }
1411  }
1412  }
1413  return selected;
1414 }
1415 
1416 //________________________________________________________________________
1417 void AliAnalysisTaskEmcal::AddObjectToEvent(TObject *obj, Bool_t attempt)
1418 {
1419  // Add object to event
1420 
1421  if (!(InputEvent()->FindListObject(obj->GetName()))) {
1422  InputEvent()->AddObject(obj);
1423  }
1424  else {
1425  if (!attempt) {
1426  AliFatal(Form("%s: Container with name %s already present. Aborting", GetName(), obj->GetName()));
1427  }
1428  }
1429 }
1430 
1431 //________________________________________________________________________
1432 Bool_t AliAnalysisTaskEmcal::IsTrackInEmcalAcceptance(AliVParticle* part, Double_t edges) const
1433 {
1434  // Determines if a track is inside the EMCal acceptance, using eta/phi at the vertex (no propagation).
1435  // Includes +/- edges. Useful to determine whether track propagation should be attempted.
1436 
1437  if (!fGeom) {
1438  AliWarning(Form("%s - AliAnalysisTaskEmcal::IsTrackInEmcalAcceptance - Geometry is not available!", GetName()));
1439  return kFALSE;
1440  }
1441 
1442  Double_t minPhi = fGeom->GetArm1PhiMin() - edges;
1443  Double_t maxPhi = fGeom->GetArm1PhiMax() + edges;
1444 
1445  if (part->Phi() > minPhi && part->Phi() < maxPhi) {
1446  return kTRUE;
1447  }
1448  else {
1449  return kFALSE;
1450  }
1451 }
1452 
1453 //________________________________________________________________________
1455 {
1456  axis->SetBinLabel(1, "NullObject");
1457  axis->SetBinLabel(2, "Pt");
1458  axis->SetBinLabel(3, "Acceptance");
1459  axis->SetBinLabel(4, "BitMap");
1460  axis->SetBinLabel(5, "Bit4");
1461  axis->SetBinLabel(6, "Bit5");
1462  axis->SetBinLabel(7, "Bit6");
1463  axis->SetBinLabel(8, "NotHybridTrack");
1464  axis->SetBinLabel(9, "MCFlag");
1465  axis->SetBinLabel(10, "MCGenerator");
1466  axis->SetBinLabel(11, "ChargeCut");
1467  axis->SetBinLabel(12, "MinDistanceTPCSectorEdge");
1468  axis->SetBinLabel(13, "MinMCLabelAccept");
1469  axis->SetBinLabel(14, "IsEMCal");
1470  axis->SetBinLabel(15, "Time");
1471  axis->SetBinLabel(16, "Energy");
1472  axis->SetBinLabel(17, "Bit16");
1473  axis->SetBinLabel(18, "Bit17");
1474  axis->SetBinLabel(19, "Area");
1475  axis->SetBinLabel(20, "AreaEmc");
1476  axis->SetBinLabel(21, "ZLeadingCh");
1477  axis->SetBinLabel(22, "ZLeadingEmc");
1478  axis->SetBinLabel(23, "NEF");
1479  axis->SetBinLabel(24, "MinLeadPt");
1480  axis->SetBinLabel(25, "MaxTrackPt");
1481  axis->SetBinLabel(26, "MaxClusterPt");
1482  axis->SetBinLabel(27, "Flavour");
1483  axis->SetBinLabel(28, "TagStatus");
1484  axis->SetBinLabel(29, "MinNConstituents");
1485  axis->SetBinLabel(30, "Bit29");
1486  axis->SetBinLabel(31, "Bit30");
1487  axis->SetBinLabel(32, "Bit31");
1488 }
1489 
1490 //________________________________________________________________________
1491 Double_t AliAnalysisTaskEmcal::GetParallelFraction(AliVParticle* part1, AliVParticle* part2)
1492 {
1493  // Calculates the fraction of momentum of part 1 w.r.t. part 2 in the direction of part 2.
1494 
1495  TVector3 vect1(part1->Px(), part1->Py(), part1->Pz());
1496  TVector3 vect2(part2->Px(), part2->Py(), part2->Pz());
1497  Double_t z = (vect1 * vect2) / (vect2 * vect2);
1498  return z;
1499 }
1500 
1501 //________________________________________________________________________
1502 Double_t AliAnalysisTaskEmcal::GetParallelFraction(const TVector3& vect1, AliVParticle* part2)
1503 {
1504  // Calculates the fraction of momentum of vect 1 w.r.t. part 2 in the direction of part 2.
1505 
1506  TVector3 vect2(part2->Px(), part2->Py(), part2->Pz());
1507  Double_t z = (vect1 * vect2) / (vect2 * vect2);
1508  return z;
1509 }
1510 
1511 //_________________________________________________________________________________________________
1512 void AliAnalysisTaskEmcal::GetEtaPhiDiff(const AliVTrack *t, const AliVCluster *v, Double_t &phidiff, Double_t &etadiff)
1513 {
1514  // Calculate phi and eta difference between track and cluster.
1515 
1516  phidiff = 999;
1517  etadiff = 999;
1518 
1519  if (!t||!v) return;
1520 
1521  Double_t veta = t->GetTrackEtaOnEMCal();
1522  Double_t vphi = t->GetTrackPhiOnEMCal();
1523 
1524  Float_t pos[3] = {0};
1525  v->GetPosition(pos);
1526  TVector3 cpos(pos);
1527  Double_t ceta = cpos.Eta();
1528  Double_t cphi = cpos.Phi();
1529  etadiff=veta-ceta;
1530  phidiff=TVector2::Phi_mpi_pi(vphi-cphi);
1531 }
1532 
1533 //_________________________________________________________________________________________________
1534 Byte_t AliAnalysisTaskEmcal::GetTrackType(const AliVTrack *t)
1535 {
1536  // Get track type encoded from bits 20 and 21.
1537 
1538  Byte_t ret = 0;
1539  if (t->TestBit(BIT(22)) && !t->TestBit(BIT(23)))
1540  ret = 1;
1541  else if (!t->TestBit(BIT(22)) && t->TestBit(BIT(23)))
1542  ret = 2;
1543  else if (t->TestBit(BIT(22)) && t->TestBit(BIT(23)))
1544  ret = 3;
1545  return ret;
1546 }
1547 
1548 //________________________________________________________________________
1549 Byte_t AliAnalysisTaskEmcal::GetTrackType(const AliAODTrack *aodTrack, UInt_t filterBit1, UInt_t filterBit2)
1550 {
1551  // Return track type: 0 = filterBit1, 1 = filterBit2 && ITS, 2 = filterBit2 && !ITS.
1552  // Returns 3 if filterBit1 and filterBit2 do not test.
1553  // WARNING: only works with AOD tracks and AOD filter bits must be provided. Otherwise will always return 0.
1554 
1555  Int_t res = 0;
1556 
1557  if (aodTrack->TestFilterBit(filterBit1)) {
1558  res = 0;
1559  }
1560  else if (aodTrack->TestFilterBit(filterBit2)) {
1561  if ((aodTrack->GetStatus()&AliVTrack::kITSrefit)!=0) {
1562  res = 1;
1563  }
1564  else {
1565  res = 2;
1566  }
1567  }
1568  else {
1569  res = 3;
1570  }
1571 
1572  return res;
1573 }
Bool_t AcceptCluster(AliVCluster *vp)
void SetParticlePtCut(Double_t cut)
TH1 * fHistTrials
x section from pythia header
Bool_t AcceptTrack(AliVParticle *track, Int_t c=0) const
Bool_t HasTriggerType(TriggerType triggersel)
Int_t fNTrials
event pt hard bin
AliEmcalTriggerPatchInfo * GetMainTriggerPatch(TriggerCategory triggersel=kTriggerLevel1Jet, Bool_t doOfflinSimple=kFALSE)
TSystem * gSystem
Double_t fPtHard
event Pythia header
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 fEPV0
event centrality bin
TList * list
virtual void SetArray(AliVEvent *event)
Bool_t AcceptCluster(AliVCluster *clus, Int_t c=0) const
Int_t fCentBin
event centrality
TH1 * fHistEventsAfterSel
total number of trials per pt hard bin after selection
Bool_t AcceptParticle(AliVParticle *vp)
void SetArrayName(const char *n)
TH1 * fHistEventPlane
z vertex position
TList * fOutput
x-section from pythia header
TH1 * fHistEvents
trials from pyxsec.root
void SetClusPtCut(Double_t cut, Int_t c=0)
AliClusterContainer * AddClusterContainer(const char *n)
Double_t fEPV0C
event plane V0A
TH1 * fHistCentrality
pt hard distribution
void SetTrackEtaLimits(Double_t min, Double_t max, Int_t c=0)
AliParticleContainer * GetParticleContainer(Int_t i=0) const
TProfile * fHistXsectionAfterSel
total number of events per pt hard bin after selection
EMCalTriggerMode_t fEMCalTriggerMode
virtual Bool_t FillHistograms()
AliVCluster * GetAcceptCluster(Int_t i)
Int_t GetNParticles(Int_t i=0) const
TClonesArray * fCaloClusters
tracks
ClassImp(AliAnalysisTaskEmcal) AliAnalysisTaskEmcal
Bool_t IsTrackInEmcalAcceptance(AliVParticle *part, Double_t edges=0.9) const
void SetArray(AliVEvent *event)
TH1 * fHistTriggerClasses
book keep reasons for rejecting event
AliEMCALGeometry * fGeom
whether it's an ESD analysis
AliVParticle * GetAcceptParticle(Int_t i)
AliGenPythiaEventHeader * fPythiaHeader
event beam type
void SetTrackPhiLimits(Double_t min, Double_t max, Int_t c=0)
AliParticleContainer * AddParticleContainer(const char *n)
AliAnalysisUtils * fAliAnalysisUtils
AliClusterContainer * GetClusterContainer(Int_t i=0) const
virtual Bool_t FillGeneralHistograms()
TClonesArray * GetParticleArray(Int_t i=0) const
BeamType fBeamType
event vertex number of contributors
Double_t fCent
trigger patch info array
TClonesArray * GetArray() const
Int_t GetNClusters(Int_t i=0) const
Int_t fNVertCont
event vertex
static Double_t GetParallelFraction(AliVParticle *part1, AliVParticle *part2)
virtual Bool_t RetrieveEventObjects()
TProfile * fHistXsection
total number of events per pt hard bin
Bool_t PythiaInfoFromFile(const char *currFile, Float_t &fXsec, Float_t &fTrials, Int_t &pthard)
void UserExec(Option_t *option)
AliVCaloCells * fCaloCells
clusters
const TString & GetArrayName() const
TClonesArray * GetArrayFromEvent(const char *name, const char *clname=0)
virtual Bool_t IsEventSelected()
TH1 * fHistPtHard
x section from pyxsec.root
void SetParticleEtaLimits(Double_t min, Double_t max)
Int_t fPtHardBin
event pt hard
TClonesArray * fTracks
emcal geometry
TH1 * fHistTrialsAfterSel
incoming and selected events
Bool_t fIsEsd
vertex selection (optional)
Double_t fVertex[3]
event plane V0C
TFile * file
TH1 * fHistEventRejection
event plane distribution
Int_t nevents[nsamples]
TClonesArray * fTriggerPatchInfo
calo triggers
TClonesArray * GetClusterArray(Int_t i=0) const
Double_t fEPV0A
event plane V0
void SetParticlePhiLimits(Double_t min, Double_t max, Double_t offset=0.)
void SetClusPtCut(Double_t cut)
void AddObjectToEvent(TObject *obj, Bool_t attempt=kFALSE)
AliVCaloTrigger * fCaloTriggers
cells
void SetRejectionReasonLabels(TAxis *axis)
TH1 * fHistZVertex
event centrality distribution
Int_t GetNEntries() const
static Byte_t GetTrackType(const AliVTrack *t)
void SetClusTimeCut(Double_t min, Double_t max, Int_t c=0)
Float_t fXsection
event trials
TH1 * fHistEventCount
output list
void SetClusTimeCut(Double_t min, Double_t max)
AliVParticle * GetAcceptParticleFromArray(Int_t p, Int_t c=0) const
AliVCluster * GetAcceptClusterFromArray(Int_t cl, Int_t c=0) const
static Double_t fgkEMCalDCalPhiDivide