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