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