AliPhysics  de71be2 (de71be2)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
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(250),
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(250),
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 particle container
1224  // will be called in AddTask macro
1225 
1226  if (TString(n).IsNull()) return 0;
1227 
1228  AliTrackContainer* cont = new AliTrackContainer(n);
1229 
1230  fParticleCollArray.Add(cont);
1231 
1232  return cont;
1233 }
1234 
1235 //________________________________________________________________________
1237 {
1238  // Add particle container
1239  // will be called in AddTask macro
1240 
1241  if (TString(n).IsNull()) return 0;
1242 
1244 
1245  fParticleCollArray.Add(cont);
1246 
1247  return cont;
1248 }
1249 
1250 //________________________________________________________________________
1252 {
1253  // Add cluster container
1254  // will be called in AddTask macro
1255 
1256  if (TString(n).IsNull()) return 0;
1257 
1259 
1260  fClusterCollArray.Add(cont);
1261 
1262  return cont;
1263 }
1264 
1265 //________________________________________________________________________
1267 {
1268  // Get i^th particle container
1269 
1270  if (i<0 || i>fParticleCollArray.GetEntriesFast()) return 0;
1271  AliParticleContainer *cont = static_cast<AliParticleContainer*>(fParticleCollArray.At(i));
1272  return cont;
1273 }
1274 
1275 //________________________________________________________________________
1277 {
1278  // Get i^th cluster container
1279 
1280  if (i<0 || i>fClusterCollArray.GetEntriesFast()) return 0;
1281  AliClusterContainer *cont = static_cast<AliClusterContainer*>(fClusterCollArray.At(i));
1282  return cont;
1283 }
1284 
1285 //________________________________________________________________________
1287 {
1288  // Get particle container with name
1289 
1290  AliParticleContainer *cont = static_cast<AliParticleContainer*>(fParticleCollArray.FindObject(name));
1291  return cont;
1292 }
1293 
1294 //________________________________________________________________________
1296 {
1297  // Get cluster container with name
1298 
1299  AliClusterContainer *cont = static_cast<AliClusterContainer*>(fClusterCollArray.FindObject(name));
1300  return cont;
1301 }
1302 
1303 //________________________________________________________________________
1304 TClonesArray* AliAnalysisTaskEmcal::GetParticleArray(Int_t i) const
1305 {
1306  // Get i^th TClonesArray with AliVParticle
1307 
1309  if (!cont) {
1310  AliError(Form("%s: Particle container %d not found",GetName(),i));
1311  return 0;
1312  }
1313  TString contName = cont->GetArrayName();
1314  return cont->GetArray();
1315 }
1316 
1317 //________________________________________________________________________
1318 TClonesArray* AliAnalysisTaskEmcal::GetClusterArray(Int_t i) const
1319 {
1320  // Get i^th TClonesArray with AliVCluster
1321 
1323  if (!cont) {
1324  AliError(Form("%s:Cluster container %d not found",GetName(),i));
1325  return 0;
1326  }
1327  return cont->GetArray();
1328 }
1329 
1330 //________________________________________________________________________
1331 AliVParticle* AliAnalysisTaskEmcal::GetAcceptParticleFromArray(Int_t p, Int_t c) const
1332 {
1333  // Get particle p if accepted from container c
1334  // If particle not accepted return 0
1335 
1337  if (!cont) {
1338  AliError(Form("%s: Particle container %d not found",GetName(),c));
1339  return 0;
1340  }
1341  AliVParticle *vp = cont->GetAcceptParticle(p);
1342 
1343  return vp;
1344 }
1345 
1346 //________________________________________________________________________
1347 AliVCluster* AliAnalysisTaskEmcal::GetAcceptClusterFromArray(Int_t cl, Int_t c) const
1348 {
1349  // Get particle p if accepted from container c
1350  // If particle not accepted return 0
1351 
1353  if (!cont) {
1354  AliError(Form("%s: Cluster container %d not found",GetName(),c));
1355  return 0;
1356  }
1357  AliVCluster *vc = cont->GetAcceptCluster(cl);
1358 
1359  return vc;
1360 }
1361 
1362 //________________________________________________________________________
1364 {
1365  // Get number of entries in particle array i
1366 
1368  if (!cont) {
1369  AliError(Form("%s: Particle container %d not found",GetName(),i));
1370  return 0;
1371  }
1372  return cont->GetNEntries();
1373 }
1374 
1375 //________________________________________________________________________
1377 {
1378  // Get number of entries in cluster array i
1379 
1381  if (!cont) {
1382  AliError(Form("%s: Cluster container %d not found",GetName(),i));
1383  return 0;
1384  }
1385  return cont->GetNEntries();
1386 }
1387 
1388 //________________________________________________________________________
1389 AliEMCALTriggerPatchInfo* AliAnalysisTaskEmcal::GetMainTriggerPatch(TriggerCategory trigger, Bool_t doSimpleOffline)
1390 {
1391  // Get main trigger match
1392  //
1393  // For the selection of the main trigger patch, high and low threshold triggers of a given category are grouped
1394  // If there are more than 1 main patch of a given trigger category (i.e. different high and low threshold patches),
1395  // the highest one according to the ADC value is taken. In case doSimpleOffline is true, then only the patches from
1396  // the simple offline trigger are used.
1397 
1398  if (!fTriggerPatchInfo) {
1399  AliError(Form("%s: fTriggerPatchInfo not available",GetName()));
1400  return 0;
1401  }
1402 
1403  //number of patches in event
1404  Int_t nPatch = fTriggerPatchInfo->GetEntries();
1405 
1406  //extract main trigger patch(es)
1407  AliEMCALTriggerPatchInfo *patch(NULL), *selected(NULL);
1408  for (Int_t iPatch = 0; iPatch < nPatch; iPatch++) {
1409 
1410  patch = (AliEMCALTriggerPatchInfo*)fTriggerPatchInfo->At( iPatch );
1411  if (patch->IsMainTrigger()) {
1412  if(doSimpleOffline){
1413  if(patch->IsOfflineSimple()){
1414  switch(trigger){
1415  case kTriggerLevel0:
1416  // option not yet implemented in the trigger maker
1417  if(patch->IsLevel0()) selected = patch;
1418  break;
1419  case kTriggerLevel1Jet:
1420  if(patch->IsJetHighSimple() || patch->IsJetLowSimple()){
1421  if(!selected) selected = patch;
1422  else if(patch->GetADCOfflineAmp() > selected->GetADCOfflineAmp()) selected = patch;
1423  }
1424  break;
1425  case kTriggerLevel1Gamma:
1426  if(patch->IsGammaHighSimple() || patch->IsGammaLowSimple()){
1427  if(!selected) selected = patch;
1428  else if(patch->GetADCOfflineAmp() > selected->GetADCOfflineAmp()) selected = patch;
1429  }
1430  break;
1431  default: // Silence compiler warnings
1432  AliError("Untreated case: Main Patch is recalculated; should be in 'else' branch");
1433  };
1434  }
1435  } else { // Not OfflineSimple
1436  switch(trigger){
1437  case kTriggerLevel0:
1438  if(patch->IsLevel0()) selected = patch;
1439  break;
1440  case kTriggerLevel1Jet:
1441  if(patch->IsJetHigh() || patch->IsJetLow()){
1442  if(!selected) selected = patch;
1443  else if (patch->GetADCAmp() > selected->GetADCAmp())
1444  selected = patch;
1445  }
1446  break;
1447  case kTriggerLevel1Gamma:
1448  if(patch->IsGammaHigh() || patch->IsGammaLow()){
1449  if(!selected) selected = patch;
1450  else if (patch->GetADCAmp() > selected->GetADCAmp())
1451  selected = patch;
1452  }
1453  break;
1454  default:
1455  AliError("Untreated case: Main Patch is recalculated; should be in 'else' branch");
1456  };
1457  }
1458  }
1459  else if ((trigger == kTriggerRecalcJet && patch->IsRecalcJet()) ||
1460  (trigger == kTriggerRecalcGamma && patch->IsRecalcGamma())) { // recalculated patches
1461  if (doSimpleOffline && patch->IsOfflineSimple()) {
1462  if(!selected) selected = patch;
1463  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
1464  selected = patch;
1465  }
1466  else if (!doSimpleOffline && !patch->IsOfflineSimple()) {
1467  if(!selected) selected = patch;
1468  else if (patch->GetADCAmp() > selected->GetADCAmp())
1469  selected = patch;
1470  }
1471  }
1472  }
1473  return selected;
1474 }
1475 
1476 //________________________________________________________________________
1477 void AliAnalysisTaskEmcal::AddObjectToEvent(TObject *obj, Bool_t attempt)
1478 {
1479  // Add object to event
1480 
1481  if (!(InputEvent()->FindListObject(obj->GetName()))) {
1482  InputEvent()->AddObject(obj);
1483  }
1484  else {
1485  if (!attempt) {
1486  AliFatal(Form("%s: Container with name %s already present. Aborting", GetName(), obj->GetName()));
1487  }
1488  }
1489 }
1490 
1491 //________________________________________________________________________
1492 Bool_t AliAnalysisTaskEmcal::IsTrackInEmcalAcceptance(AliVParticle* part, Double_t edges) const
1493 {
1494  // Determines if a track is inside the EMCal acceptance, using eta/phi at the vertex (no propagation).
1495  // Includes +/- edges. Useful to determine whether track propagation should be attempted.
1496 
1497  if (!fGeom) {
1498  AliWarning(Form("%s - AliAnalysisTaskEmcal::IsTrackInEmcalAcceptance - Geometry is not available!", GetName()));
1499  return kFALSE;
1500  }
1501 
1502  Double_t minPhi = fGeom->GetArm1PhiMin() - edges;
1503  Double_t maxPhi = fGeom->GetArm1PhiMax() + edges;
1504 
1505  if (part->Phi() > minPhi && part->Phi() < maxPhi) {
1506  return kTRUE;
1507  }
1508  else {
1509  return kFALSE;
1510  }
1511 }
1512 
1513 //________________________________________________________________________
1515 {
1516  axis->SetBinLabel(1, "NullObject");
1517  axis->SetBinLabel(2, "Pt");
1518  axis->SetBinLabel(3, "Acceptance");
1519  axis->SetBinLabel(4, "MCLabel");
1520  axis->SetBinLabel(5, "BitMap");
1521  axis->SetBinLabel(6, "HF cut");
1522  axis->SetBinLabel(7, "Bit6");
1523  axis->SetBinLabel(8, "NotHybridTrack");
1524  axis->SetBinLabel(9, "MCFlag");
1525  axis->SetBinLabel(10, "MCGenerator");
1526  axis->SetBinLabel(11, "ChargeCut");
1527  axis->SetBinLabel(12, "MinDistanceTPCSectorEdge");
1528  axis->SetBinLabel(13, "Bit12");
1529  axis->SetBinLabel(14, "IsEMCal");
1530  axis->SetBinLabel(15, "Time");
1531  axis->SetBinLabel(16, "Energy");
1532  axis->SetBinLabel(17, "ExoticCut");
1533  axis->SetBinLabel(18, "Bit17");
1534  axis->SetBinLabel(19, "Area");
1535  axis->SetBinLabel(20, "AreaEmc");
1536  axis->SetBinLabel(21, "ZLeadingCh");
1537  axis->SetBinLabel(22, "ZLeadingEmc");
1538  axis->SetBinLabel(23, "NEF");
1539  axis->SetBinLabel(24, "MinLeadPt");
1540  axis->SetBinLabel(25, "MaxTrackPt");
1541  axis->SetBinLabel(26, "MaxClusterPt");
1542  axis->SetBinLabel(27, "Flavour");
1543  axis->SetBinLabel(28, "TagStatus");
1544  axis->SetBinLabel(29, "MinNConstituents");
1545  axis->SetBinLabel(30, "Bit29");
1546  axis->SetBinLabel(31, "Bit30");
1547  axis->SetBinLabel(32, "Bit31");
1548 }
1549 
1550 //________________________________________________________________________
1551 Double_t AliAnalysisTaskEmcal::GetParallelFraction(AliVParticle* part1, AliVParticle* part2)
1552 {
1553  // Calculates the fraction of momentum of part 1 w.r.t. part 2 in the direction of part 2.
1554 
1555  TVector3 vect1(part1->Px(), part1->Py(), part1->Pz());
1556  TVector3 vect2(part2->Px(), part2->Py(), part2->Pz());
1557  Double_t z = (vect1 * vect2) / (vect2 * vect2);
1558  return z;
1559 }
1560 
1561 //________________________________________________________________________
1562 Double_t AliAnalysisTaskEmcal::GetParallelFraction(const TVector3& vect1, AliVParticle* part2)
1563 {
1564  // Calculates the fraction of momentum of vect 1 w.r.t. part 2 in the direction of part 2.
1565 
1566  TVector3 vect2(part2->Px(), part2->Py(), part2->Pz());
1567  Double_t z = (vect1 * vect2) / (vect2 * vect2);
1568  return z;
1569 }
1570 
1571 //_________________________________________________________________________________________________
1572 void AliAnalysisTaskEmcal::GetEtaPhiDiff(const AliVTrack *t, const AliVCluster *v, Double_t &phidiff, Double_t &etadiff)
1573 {
1574  // Calculate phi and eta difference between track and cluster.
1575 
1576  phidiff = 999;
1577  etadiff = 999;
1578 
1579  if (!t||!v) return;
1580 
1581  Double_t veta = t->GetTrackEtaOnEMCal();
1582  Double_t vphi = t->GetTrackPhiOnEMCal();
1583 
1584  Float_t pos[3] = {0};
1585  v->GetPosition(pos);
1586  TVector3 cpos(pos);
1587  Double_t ceta = cpos.Eta();
1588  Double_t cphi = cpos.Phi();
1589  etadiff=veta-ceta;
1590  phidiff=TVector2::Phi_mpi_pi(vphi-cphi);
1591 }
1592 
1593 //_________________________________________________________________________________________________
1594 Byte_t AliAnalysisTaskEmcal::GetTrackType(const AliVTrack *t)
1595 {
1596  // Get track type encoded from bits 20 and 21.
1597 
1598  Byte_t ret = 0;
1599  if (t->TestBit(BIT(22)) && !t->TestBit(BIT(23)))
1600  ret = 1;
1601  else if (!t->TestBit(BIT(22)) && t->TestBit(BIT(23)))
1602  ret = 2;
1603  else if (t->TestBit(BIT(22)) && t->TestBit(BIT(23)))
1604  ret = 3;
1605  return ret;
1606 }
1607 
1608 //________________________________________________________________________
1609 Byte_t AliAnalysisTaskEmcal::GetTrackType(const AliAODTrack *aodTrack, UInt_t filterBit1, UInt_t filterBit2)
1610 {
1611  // Return track type: 0 = filterBit1, 1 = filterBit2 && ITS, 2 = filterBit2 && !ITS.
1612  // Returns 3 if filterBit1 and filterBit2 do not test.
1613  // WARNING: only works with AOD tracks and AOD filter bits must be provided. Otherwise will always return 0.
1614 
1615  Int_t res = 0;
1616 
1617  if (aodTrack->TestFilterBit(filterBit1)) {
1618  res = 0;
1619  }
1620  else if (aodTrack->TestFilterBit(filterBit2)) {
1621  if ((aodTrack->GetStatus()&AliVTrack::kITSrefit)!=0) {
1622  res = 1;
1623  }
1624  else {
1625  res = 2;
1626  }
1627  }
1628  else {
1629  res = 3;
1630  }
1631 
1632  return res;
1633 }
1634 
1635 //________________________________________________________________________
1637 {
1638  // Copy some information about the Pythia event in a PythaInfo object
1639 
1640  if (!fPythiaInfo) {
1642  }
1643 
1644  AliStack* stack = mcEvent->Stack();
1645 
1646  const Int_t nprim = stack->GetNprimary();
1647  // reject if partons are missing from stack for some reason
1648  if (nprim < 8) return;
1649 
1650  TParticle *part6 = stack->Particle(6);
1651  TParticle *part7 = stack->Particle(7);
1652 
1653  fPythiaInfo->SetPartonFlag6(TMath::Abs(part6->GetPdgCode()));
1654  fPythiaInfo->SetParton6(part6->Pt(), part6->Eta(), part6->Phi(), part6->GetMass());
1655 
1656  fPythiaInfo->SetPartonFlag7(TMath::Abs(part7->GetPdgCode()));
1657  fPythiaInfo->SetParton7(part7->Pt(), part7->Eta(), part7->Phi(), part7->GetMass());
1658 
1659  AliGenPythiaEventHeader *pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(mcEvent->GenEventHeader());
1660  if(pythiaGenHeader){
1661  Float_t ptWeight=pythiaGenHeader->EventWeight();
1662  fPythiaInfo->SetPythiaEventWeight(ptWeight);}
1663 
1664 
1665 }
void SetParticlePtCut(Double_t cut)
void SetParton7(Float_t pt, Float_t eta, Float_t phi, Float_t mass=0)
TH1 * fHistTrials
x section from pythia header
AliEmcalPythiaInfo * fPythiaInfo
x-section from pythia header
Bool_t AcceptTrack(AliVParticle *track, Int_t c=0) const
Bool_t HasTriggerType(TriggerType triggersel)
Int_t fNTrials
event pt hard bin
void SetPartonFlag7(Int_t flag7)
TSystem * gSystem
Double_t fPtHard
event Pythia header
void SetTrackPtCut(Double_t cut, Int_t c=0)
static void GetEtaPhiDiff(const AliVTrack *t, const AliVCluster *v, Double_t &phidiff, Double_t &etadiff)
Double_t fEPV0
event centrality bin
TList * list
virtual void SetArray(AliVEvent *event)
Bool_t AcceptCluster(AliVCluster *clus, Int_t c=0) const
Int_t fCentBin
event centrality
TH1 * fHistEventsAfterSel
total number of trials per pt hard bin after selection
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
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
AliMCParticleContainer * AddMCParticleContainer(const char *n)
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
AliTrackContainer * AddTrackContainer(const char *n)
Store some informaion about a Pythia eventThis class is used to store some information about a Pythia...
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
virtual Bool_t AcceptCluster(Int_t i)
Declaration of class AliEmcalPythiaInfo.
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
virtual AliVParticle * GetAcceptParticle(Int_t i=-1)
void SetClusTimeCut(Double_t min, Double_t max)
void SetParticlePhiLimits(Double_t min, Double_t max)
AliVParticle * GetAcceptParticleFromArray(Int_t p, Int_t c=0) const
virtual Bool_t AcceptParticle(const AliVParticle *vp)
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