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