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