]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWG/EMCAL/AliAnalysisTaskEmcal.cxx
Add ptr to original AliVTrack.
[u/mrichter/AliRoot.git] / PWG / EMCAL / AliAnalysisTaskEmcal.cxx
1 // $Id$
2 //
3 // Emcal base analysis task.
4 //
5 // Author: S.Aiola, M. Verweij
6
7 #include "AliAnalysisTaskEmcal.h"
8
9 #include <TClonesArray.h>
10 #include <TList.h>
11 #include <TObject.h>
12 #include <TH1F.h>
13 #include <TProfile.h>
14 #include <TSystem.h>
15 #include <TFile.h>
16 #include <TChain.h>
17 #include <TKey.h>
18
19 #include "AliAODEvent.h"
20 #include "AliAnalysisManager.h"
21 #include "AliCentrality.h"
22 #include "AliEMCALGeometry.h"
23 #include "AliESDEvent.h"
24 #include "AliEmcalParticle.h"
25 #include "AliEventplane.h"
26 #include "AliInputEventHandler.h"
27 #include "AliLog.h"
28 #include "AliMCParticle.h"
29 #include "AliVCluster.h"
30 #include "AliVEventHandler.h"
31 #include "AliVParticle.h"
32 #include "AliVCaloTrigger.h"
33 #include "AliGenPythiaEventHeader.h"
34 #include "AliAODMCHeader.h"
35 #include "AliMCEvent.h"
36 #include "AliAnalysisUtils.h"
37 #include "AliEmcalTriggerPatchInfo.h"
38
39 #include "AliParticleContainer.h"
40 #include "AliClusterContainer.h"
41
42 ClassImp(AliAnalysisTaskEmcal)
43
44 //________________________________________________________________________
45 AliAnalysisTaskEmcal::AliAnalysisTaskEmcal() : 
46   AliAnalysisTaskSE("AliAnalysisTaskEmcal"),
47   fForceBeamType(kNA),
48   fGeneralHistograms(kFALSE),
49   fInitialized(kFALSE),
50   fCreateHisto(kTRUE),
51   fCaloCellsName(),
52   fCaloTriggersName(),
53   fCaloTriggerPatchInfoName(),
54   fMinCent(-999),
55   fMaxCent(-999),
56   fMinVz(-999),
57   fMaxVz(-999),
58   fTrackPtCut(0),
59   fMinNTrack(0),
60   fUseAliAnaUtils(kFALSE),
61   fAliAnalysisUtils(0x0),
62   fOffTrigger(AliVEvent::kAny),
63   fTrigClass(),
64   fTriggerTypeSel(kND),
65   fNbins(500),
66   fMinBinPt(0),
67   fMaxBinPt(250),
68   fMinPtTrackInEmcal(0),
69   fEventPlaneVsEmcal(-1),
70   fMinEventPlane(-1e6),
71   fMaxEventPlane(1e6),
72   fCentEst("V0M"),
73   fIsEmbedded(kFALSE),
74   fIsPythia(kFALSE),
75   fSelectPtHardBin(-999),
76   fMinMCLabel(0),
77   fMCLabelShift(0),
78   fNcentBins(4),
79   fGeom(0),
80   fTracks(0),
81   fCaloClusters(0),
82   fCaloCells(0),
83   fCaloTriggers(0),
84   fTriggerPatchInfo(0),
85   fCent(0),
86   fCentBin(-1),
87   fEPV0(-1.0),
88   fEPV0A(-1.0),
89   fEPV0C(-1.0),
90   fNVertCont(0),
91   fBeamType(kNA),
92   fPythiaHeader(0),
93   fPtHard(0),
94   fPtHardBin(0),
95   fNTrials(0),
96   fParticleCollArray(),
97   fClusterCollArray(),
98   fMainTriggerPatch(0x0),
99   fTriggerType(kND),
100   fOutput(0),
101   fHistTrialsAfterSel(0),
102   fHistEventsAfterSel(0),
103   fHistTrials(0),
104   fHistXsection(0),
105   fHistEvents(0),
106   fHistPtHard(0),
107   fHistCentrality(0),
108   fHistZVertex(0),
109   fHistEventPlane(0),
110   fHistEventRejection(0)
111 {
112   // Default constructor.
113
114   fVertex[0] = 0;
115   fVertex[1] = 0;
116   fVertex[2] = 0;
117
118   fParticleCollArray.SetOwner(kTRUE);
119   fClusterCollArray.SetOwner(kTRUE);
120 }
121
122 //________________________________________________________________________
123 AliAnalysisTaskEmcal::AliAnalysisTaskEmcal(const char *name, Bool_t histo) : 
124   AliAnalysisTaskSE(name),
125   fForceBeamType(kNA),
126   fGeneralHistograms(kFALSE),
127   fInitialized(kFALSE),
128   fCreateHisto(histo),
129   fCaloCellsName(),
130   fCaloTriggersName(),
131   fCaloTriggerPatchInfoName(),
132   fMinCent(-999),
133   fMaxCent(-999),
134   fMinVz(-999),
135   fMaxVz(-999),
136   fTrackPtCut(0),
137   fMinNTrack(0),
138   fUseAliAnaUtils(kFALSE),
139   fAliAnalysisUtils(0x0),
140   fOffTrigger(AliVEvent::kAny),
141   fTrigClass(),
142   fTriggerTypeSel(kND),
143   fNbins(500),
144   fMinBinPt(0),
145   fMaxBinPt(250),
146   fMinPtTrackInEmcal(0),
147   fEventPlaneVsEmcal(-1),
148   fMinEventPlane(-1e6),
149   fMaxEventPlane(1e6),
150   fCentEst("V0M"),
151   fIsEmbedded(kFALSE),
152   fIsPythia(kFALSE),
153   fSelectPtHardBin(-999),
154   fMinMCLabel(0),
155   fMCLabelShift(0),
156   fNcentBins(4),
157   fGeom(0),
158   fTracks(0),
159   fCaloClusters(0),
160   fCaloCells(0),
161   fCaloTriggers(0),
162   fTriggerPatchInfo(0),
163   fCent(0),
164   fCentBin(-1),
165   fEPV0(-1.0),
166   fEPV0A(-1.0),
167   fEPV0C(-1.0),
168   fNVertCont(0),
169   fBeamType(kNA),
170   fPythiaHeader(0),
171   fPtHard(0),
172   fPtHardBin(0),
173   fNTrials(0),
174   fParticleCollArray(),
175   fClusterCollArray(),
176   fMainTriggerPatch(0x0),
177   fTriggerType(kND),
178   fOutput(0),
179   fHistTrialsAfterSel(0),
180   fHistEventsAfterSel(0),
181   fHistTrials(0),
182   fHistXsection(0),
183   fHistEvents(0),
184   fHistPtHard(0),
185   fHistCentrality(0),
186   fHistZVertex(0),
187   fHistEventPlane(0),
188   fHistEventRejection(0)
189 {
190   // Standard constructor.
191
192   fVertex[0] = 0;
193   fVertex[1] = 0;
194   fVertex[2] = 0;
195
196   fParticleCollArray.SetOwner(kTRUE);
197   fClusterCollArray.SetOwner(kTRUE);
198
199   if (fCreateHisto) {
200     DefineOutput(1, TList::Class()); 
201   }
202 }
203
204 //________________________________________________________________________
205 AliAnalysisTaskEmcal::~AliAnalysisTaskEmcal()
206 {
207   // Destructor
208 }
209
210 //________________________________________________________________________
211 void AliAnalysisTaskEmcal::SetClusPtCut(Double_t cut, Int_t c)
212 {
213   AliClusterContainer *cont = GetClusterContainer(c);
214   if (cont) cont->SetClusPtCut(cut);
215   else AliError(Form("%s in SetClusPtCut(...): container %d not found",GetName(),c));
216 }
217
218 //________________________________________________________________________
219 void AliAnalysisTaskEmcal::SetClusTimeCut(Double_t min, Double_t max, Int_t c)
220 {
221   AliClusterContainer *cont = GetClusterContainer(c);
222   if (cont) cont->SetClusTimeCut(min,max);
223   else AliError(Form("%s in SetClusTimeCut(...): container %d not found",GetName(),c));
224 }
225
226 //________________________________________________________________________
227 void AliAnalysisTaskEmcal::SetTrackPtCut(Double_t cut, Int_t c)
228 {
229   AliParticleContainer *cont = GetParticleContainer(c);
230   if (cont) cont->SetParticlePtCut(cut);
231   else AliError(Form("%s in SetTrackPtCut(...): container %d not found",GetName(),c));
232
233   fTrackPtCut = cut;
234 }
235
236 //________________________________________________________________________
237 void AliAnalysisTaskEmcal::SetTrackEtaLimits(Double_t min, Double_t max, Int_t c)
238 {
239   AliParticleContainer *cont = GetParticleContainer(c);
240   if (cont) cont->SetParticleEtaLimits(min,max);
241   else AliError(Form("%s in SetTrackPtCut(...): container %d not found",GetName(),c));
242 }
243
244 //________________________________________________________________________
245 void AliAnalysisTaskEmcal::SetTrackPhiLimits(Double_t min, Double_t max, Int_t c)
246 {
247   AliParticleContainer *cont = GetParticleContainer(c);
248   if (cont) cont->SetParticlePhiLimits(min,max);
249   else AliError(Form("%s in SetTrackPhiLimits(...): container %d not found",GetName(),c));
250 }
251
252 //________________________________________________________________________
253 void AliAnalysisTaskEmcal::UserCreateOutputObjects()
254 {
255   // Create user output.
256   if (!fCreateHisto)
257     return;
258
259   OpenFile(1);
260   fOutput = new TList();
261   fOutput->SetOwner();
262
263   if (fForceBeamType == kpp)
264     fNcentBins = 1;
265
266   if (!fGeneralHistograms)
267     return;
268
269   if (fIsPythia) {
270     fHistTrialsAfterSel = new TH1F("fHistTrialsAfterSel", "fHistTrialsAfterSel", 11, 0, 11);
271     fHistTrialsAfterSel->GetXaxis()->SetTitle("p_{T} hard bin");
272     fHistTrialsAfterSel->GetYaxis()->SetTitle("trials");
273     fOutput->Add(fHistTrialsAfterSel);
274     
275     fHistEventsAfterSel = new TH1F("fHistEventsAfterSel", "fHistEventsAfterSel", 11, 0, 11);
276     fHistEventsAfterSel->GetXaxis()->SetTitle("p_{T} hard bin");
277     fHistEventsAfterSel->GetYaxis()->SetTitle("total events");
278     fOutput->Add(fHistEventsAfterSel);
279     
280     fHistTrials = new TH1F("fHistTrials", "fHistTrials", 11, 0, 11);
281     fHistTrials->GetXaxis()->SetTitle("p_{T} hard bin");
282     fHistTrials->GetYaxis()->SetTitle("trials");
283     fOutput->Add(fHistTrials);
284
285     fHistXsection = new TProfile("fHistXsection", "fHistXsection", 11, 0, 11);
286     fHistXsection->GetXaxis()->SetTitle("p_{T} hard bin");
287     fHistXsection->GetYaxis()->SetTitle("xsection");
288     fOutput->Add(fHistXsection);
289
290     fHistEvents = new TH1F("fHistEvents", "fHistEvents", 11, 0, 11);
291     fHistEvents->GetXaxis()->SetTitle("p_{T} hard bin");
292     fHistEvents->GetYaxis()->SetTitle("total events");
293     fOutput->Add(fHistEvents);
294
295     const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
296     const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
297     
298     for (Int_t i = 1; i < 12; i++) {
299       fHistTrialsAfterSel->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
300       fHistEventsAfterSel->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
301       
302       fHistTrials->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
303       fHistXsection->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
304       fHistEvents->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
305     }
306
307     fHistPtHard = new TH1F("fHistPtHard", "fHistPtHard", fNbins*2, fMinBinPt, fMaxBinPt*4);
308     fHistPtHard->GetXaxis()->SetTitle("p_{T,hard} (GeV/c)");
309     fHistPtHard->GetYaxis()->SetTitle("counts");
310     fOutput->Add(fHistPtHard);
311   }
312
313   fHistCentrality = new TH1F("fHistCentrality","Event centrality distribution", 200, 0, 100);
314   fHistCentrality->GetXaxis()->SetTitle("Centrality (%)");
315   fHistCentrality->GetYaxis()->SetTitle("counts");
316   fOutput->Add(fHistCentrality);
317
318   fHistZVertex = new TH1F("fHistZVertex","Z vertex position", 60, -30, 30);
319   fHistZVertex->GetXaxis()->SetTitle("z");
320   fHistZVertex->GetYaxis()->SetTitle("counts");
321   fOutput->Add(fHistZVertex);
322
323   fHistEventPlane = new TH1F("fHistEventPlane","Event plane", 120, -TMath::Pi(), TMath::Pi());
324   fHistEventPlane->GetXaxis()->SetTitle("event plane");
325   fHistEventPlane->GetYaxis()->SetTitle("counts");
326   fOutput->Add(fHistEventPlane);
327
328   fHistEventRejection = new TH1I("fHistEventRejection","Reasons to reject event",20,0,20);
329   fHistEventRejection->Fill("PhysSel",0);
330   fHistEventRejection->Fill("trigger",0);
331   fHistEventRejection->Fill("trigTypeSel",0);
332   fHistEventRejection->Fill("Cent",0);
333   fHistEventRejection->Fill("vertex contr.",0);
334   fHistEventRejection->Fill("Vz",0);
335   fHistEventRejection->Fill("trackInEmcal",0);
336   fHistEventRejection->Fill("minNTrack",0);
337   fHistEventRejection->Fill("VtxSel2013pA",0);
338   fHistEventRejection->Fill("PileUp",0);
339   fHistEventRejection->Fill("EvtPlane",0);
340   fHistEventRejection->Fill("SelPtHardBin",0);
341   fHistEventRejection->GetYaxis()->SetTitle("counts");
342   fOutput->Add(fHistEventRejection);
343
344   PostData(1, fOutput);
345 }
346
347 //________________________________________________________________________
348 Bool_t AliAnalysisTaskEmcal::FillGeneralHistograms()
349 {
350   if (fIsPythia) {
351     fHistEventsAfterSel->SetBinContent(fPtHardBin + 1, fHistEventsAfterSel->GetBinContent(fPtHardBin + 1) + 1);
352     fHistTrialsAfterSel->SetBinContent(fPtHardBin + 1, fHistTrialsAfterSel->GetBinContent(fPtHardBin + 1) + fNTrials);
353     fHistPtHard->Fill(fPtHard);
354   }
355
356   fHistCentrality->Fill(fCent);
357   fHistZVertex->Fill(fVertex[2]);
358   fHistEventPlane->Fill(fEPV0);
359
360   return kTRUE;
361 }
362
363 //________________________________________________________________________
364 void AliAnalysisTaskEmcal::UserExec(Option_t *) 
365 {
366   // Main loop, called for each event.
367
368   fMainTriggerPatch = NULL;
369
370   if (!fInitialized)
371     ExecOnce();
372
373   if (!fInitialized)
374     return;
375
376   if (!RetrieveEventObjects())
377     return;
378
379   if (!IsEventSelected()) 
380     return;
381
382   if (fGeneralHistograms && fCreateHisto) {
383     if (!FillGeneralHistograms())
384       return;
385   }
386
387   if (!Run())
388     return;
389
390   if (fCreateHisto) {
391     if (!FillHistograms())
392       return;
393   }
394     
395   if (fCreateHisto && fOutput) {
396     // information for this iteration of the UserExec in the container
397     PostData(1, fOutput);
398   }
399 }
400
401 //________________________________________________________________________
402 Bool_t AliAnalysisTaskEmcal::AcceptCluster(AliVCluster *clus, Int_t c) const
403 {
404   // Return true if cluster is accepted.
405
406   if (!clus)
407     return kFALSE;
408
409   AliClusterContainer *cont = GetClusterContainer(c);
410   if (!cont) {
411     AliError(Form("%s:Container %d not found",GetName(),c));
412     return 0;
413   }
414
415   return cont->AcceptCluster(clus);
416
417 }
418
419 //________________________________________________________________________
420 Bool_t AliAnalysisTaskEmcal::AcceptTrack(AliVParticle *track, Int_t c) const
421 {
422   // Return true if track is accepted.
423
424   if (!track)
425     return kFALSE;
426
427   AliParticleContainer *cont = GetParticleContainer(c);
428   if (!cont) {
429     AliError(Form("%s:Container %d not found",GetName(),c));
430     return 0;
431   }
432
433   return cont->AcceptParticle(track);
434 }
435
436 //________________________________________________________________________
437 Bool_t AliAnalysisTaskEmcal::PythiaInfoFromFile(const char* currFile, Float_t &fXsec, Float_t &fTrials, Int_t &pthard)
438 {
439   //
440   // Get the cross section and the trails either from pyxsec.root or from pysec_hists.root
441   // Get the pt hard bin from the file path
442   // This is to called in Notify and should provide the path to the AOD/ESD file
443   // (Partially copied from AliAnalysisHelperJetTasks)
444
445   TString file(currFile);  
446   fXsec = 0;
447   fTrials = 1;
448
449   if (file.Contains(".zip#")) {
450     Ssiz_t pos1 = file.Index("root_archive",12,0,TString::kExact);
451     Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
452     Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
453     file.Replace(pos+1,pos2-pos1,"");
454   }
455   else {
456     // not an archive take the basename....
457     file.ReplaceAll(gSystem->BaseName(file.Data()),"");
458   }
459   AliDebug(1,Form("File name: %s",file.Data()));
460
461   // Get the pt hard bin
462   TString strPthard(file);
463
464   strPthard.Remove(strPthard.Last('/'));
465   strPthard.Remove(strPthard.Last('/'));
466   if (strPthard.Contains("AOD")) strPthard.Remove(strPthard.Last('/'));    
467   strPthard.Remove(0,strPthard.Last('/')+1);
468   if (strPthard.IsDec()) 
469     pthard = strPthard.Atoi();
470   else 
471     AliWarning(Form("Could not extract file number from path %s", strPthard.Data()));
472
473   // 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
474   TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")); 
475   
476   if (!fxsec) {
477     // next trial fetch the histgram file
478     fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
479     if (!fxsec) {
480         // not a severe condition but inciate that we have no information
481       return kFALSE;
482     }
483     else{
484       // find the tlist we want to be independtent of the name so use the Tkey
485       TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0); 
486       if (!key) {
487         fxsec->Close();
488         return kFALSE;
489       }
490       TList *list = dynamic_cast<TList*>(key->ReadObj());
491       if (!list) {
492         fxsec->Close();
493         return kFALSE;
494       }
495       fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
496       fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
497       fxsec->Close();
498     }
499   } // no tree pyxsec.root
500   else {
501     TTree *xtree = (TTree*)fxsec->Get("Xsection");
502     if (!xtree) {
503       fxsec->Close();
504       return kFALSE;
505     }
506     UInt_t   ntrials  = 0;
507     Double_t  xsection  = 0;
508     xtree->SetBranchAddress("xsection",&xsection);
509     xtree->SetBranchAddress("ntrials",&ntrials);
510     xtree->GetEntry(0);
511     fTrials = ntrials;
512     fXsec = xsection;
513     fxsec->Close();
514   }
515   return kTRUE;
516 }
517
518 //________________________________________________________________________
519 Bool_t AliAnalysisTaskEmcal::UserNotify()
520 {
521   if (!fIsPythia || !fGeneralHistograms || !fCreateHisto)
522     return kTRUE;
523
524   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
525   if (!tree) {
526     AliError(Form("%s - UserNotify: No current tree!",GetName()));
527     return kFALSE;
528   }
529
530   Float_t xsection = 0;
531   Float_t trials   = 0;
532   Int_t   pthard   = 0;
533
534   TFile *curfile = tree->GetCurrentFile();
535   if (!curfile) {
536     AliError(Form("%s - UserNotify: No current file!",GetName()));
537     return kFALSE;
538   }
539
540   TChain *chain = dynamic_cast<TChain*>(tree);
541   if (chain)
542     tree = chain->GetTree();
543
544   Int_t nevents = tree->GetEntriesFast();
545
546   PythiaInfoFromFile(curfile->GetName(), xsection, trials, pthard);
547
548   // TODO: Workaround
549   if ((pthard < 0) || (pthard > 10))
550     pthard = 0;
551   fHistTrials->Fill(pthard, trials);
552   fHistXsection->Fill(pthard, xsection);
553   fHistEvents->Fill(pthard, nevents);
554
555   return kTRUE;
556 }
557
558 //________________________________________________________________________
559 void AliAnalysisTaskEmcal::ExecOnce()
560 {
561   // Init the analysis.
562
563   if (!InputEvent()) {
564     AliError(Form("%s: Could not retrieve event! Returning!", GetName()));
565     return;
566   }
567
568   fGeom = AliEMCALGeometry::GetInstance();
569   if (!fGeom) {
570     AliError(Form("%s: Can not create geometry", GetName()));
571     return;
572   }
573
574   if (fEventPlaneVsEmcal >= 0) {
575     Double_t ep = (fGeom->GetArm1PhiMax() + fGeom->GetArm1PhiMin()) / 2 * TMath::DegToRad() + fEventPlaneVsEmcal - TMath::Pi();
576     fMinEventPlane = ep - TMath::Pi() / 4;
577     fMaxEventPlane = ep + TMath::Pi() / 4;
578   }
579
580   //Load all requested track branches - each container knows name already
581   for (Int_t i =0; i<fParticleCollArray.GetEntriesFast(); i++) {
582     AliParticleContainer *cont = static_cast<AliParticleContainer*>(fParticleCollArray.At(i));
583     cont->SetArray(InputEvent());
584   }
585
586   if (fParticleCollArray.GetEntriesFast()>0) {
587     fTracks = GetParticleArray(0);
588     if (!fTracks) {
589       AliError(Form("%s: Could not retrieve first track branch!", GetName()));
590       return;
591     }
592   }
593
594   //Load all requested cluster branches - each container knows name already
595   for (Int_t i =0; i<fClusterCollArray.GetEntriesFast(); i++) {
596     AliClusterContainer *cont = static_cast<AliClusterContainer*>(fClusterCollArray.At(i));
597     cont->SetArray(InputEvent());
598   }
599
600   if (fClusterCollArray.GetEntriesFast()>0) {
601     fCaloClusters = GetClusterArray(0);
602     if (!fCaloClusters) {
603       AliError(Form("%s: Could not retrieve first cluster branch!", GetName()));
604       return;
605     }
606   }
607
608   if (!fCaloCellsName.IsNull() && !fCaloCells) {
609     fCaloCells =  dynamic_cast<AliVCaloCells*>(InputEvent()->FindListObject(fCaloCellsName));
610     if (!fCaloCells) {
611       AliError(Form("%s: Could not retrieve cells %s!", GetName(), fCaloCellsName.Data())); 
612       return;
613     }
614   }
615
616   if (!fCaloTriggersName.IsNull() && !fCaloTriggers) {
617     fCaloTriggers =  dynamic_cast<AliVCaloTrigger*>(InputEvent()->FindListObject(fCaloTriggersName));
618     if (!fCaloTriggers) {
619       AliError(Form("%s: Could not retrieve calo triggers %s!", GetName(), fCaloTriggersName.Data())); 
620       return;
621     }
622   }
623
624   if (!fCaloTriggerPatchInfoName.IsNull() && !fTriggerPatchInfo) {
625     fTriggerPatchInfo = GetArrayFromEvent(fCaloTriggerPatchInfoName.Data(),"AliEmcalTriggerPatchInfo");
626     if (!fTriggerPatchInfo) {
627       AliError(Form("%s: Could not retrieve calo trigger patch info %s!", GetName(), fCaloTriggerPatchInfoName.Data())); 
628       return;
629     }
630
631   }
632
633   fInitialized = kTRUE;
634 }
635
636 //_____________________________________________________
637 AliAnalysisTaskEmcal::BeamType AliAnalysisTaskEmcal::GetBeamType()
638 {
639   // Get beam type : pp-AA-pA
640   // ESDs have it directly, AODs get it from hardcoded run number ranges
641
642   if (fForceBeamType != kNA)
643     return fForceBeamType;
644
645   AliESDEvent *esd = dynamic_cast<AliESDEvent*>(InputEvent());
646   if (esd) {
647     const AliESDRun *run = esd->GetESDRun();
648     TString beamType = run->GetBeamType();
649     if (beamType == "p-p")
650       return kpp;
651     else if (beamType == "A-A")
652       return kAA;
653     else if (beamType == "p-A")
654       return kpA;
655     else
656       return kNA;
657   } else {
658     Int_t runNumber = InputEvent()->GetRunNumber();
659     if ((runNumber >= 136851 && runNumber <= 139517) ||  // LHC10h
660         (runNumber >= 166529 && runNumber <= 170593))    // LHC11h
661     { 
662       return kAA;
663     } 
664     else if ((runNumber>=188365 && runNumber <= 188366) || // LHC12g
665              (runNumber >= 195344 && runNumber <= 196608))  // LHC13b-f
666     {
667       return kpA;
668     } else {
669       return kpp;
670     }
671   }  
672 }
673
674 //_____________________________________________________
675 AliAnalysisTaskEmcal::TriggerType AliAnalysisTaskEmcal::GetTriggerType()
676 {
677   // Get trigger type: kND, kJ1, kJ2
678  
679   if (!fTriggerPatchInfo)
680     return kND;
681   
682   //number of patches in event
683   Int_t nPatch = fTriggerPatchInfo->GetEntries();
684
685   //loop over patches to define trigger type of event
686   Int_t nJ1 = 0;
687   Int_t nJ2 = 0;
688   AliEmcalTriggerPatchInfo *patch;
689   for (Int_t iPatch = 0; iPatch < nPatch; iPatch++) {
690     patch = (AliEmcalTriggerPatchInfo*)fTriggerPatchInfo->At( iPatch );
691     if (patch->IsJetHigh()) nJ1++;
692     if (patch->IsJetLow())  nJ2++;
693   }
694
695   if (nJ1>0) 
696     return kJ1;
697   else if (nJ2>0)
698     return kJ2;
699   else
700     return kND;
701  
702 }
703
704 //________________________________________________________________________
705 Bool_t AliAnalysisTaskEmcal::IsEventSelected()
706 {
707   // Check if event is selected
708
709   if (fOffTrigger != AliVEvent::kAny) {
710     UInt_t res = 0;
711     const AliESDEvent *eev = dynamic_cast<const AliESDEvent*>(InputEvent());
712     if (eev) {
713       res = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
714     } else {
715       const AliAODEvent *aev = dynamic_cast<const AliAODEvent*>(InputEvent());
716       if (aev) {
717         res = aev->GetHeader()->GetOfflineTrigger();
718       }
719     }
720     if ((res & fOffTrigger) == 0) {
721       if (fGeneralHistograms) 
722         fHistEventRejection->Fill("PhysSel",1);
723       return kFALSE;
724     }
725   }
726
727   if (!fTrigClass.IsNull()) {
728     TString fired;
729     const AliESDEvent *eev = dynamic_cast<const AliESDEvent*>(InputEvent());
730     if (eev) {
731       fired = eev->GetFiredTriggerClasses();
732     } else {
733       const AliAODEvent *aev = dynamic_cast<const AliAODEvent*>(InputEvent());
734       if (aev) {
735         fired = aev->GetFiredTriggerClasses();
736       }
737     }
738     if (!fired.Contains("-B-")) {
739       if (fGeneralHistograms) 
740         fHistEventRejection->Fill("trigger",1);
741       return kFALSE;
742     }
743     TObjArray *arr = fTrigClass.Tokenize("|");
744     if (!arr) {
745       if (fGeneralHistograms) 
746         fHistEventRejection->Fill("trigger",1);
747       return kFALSE;
748     }
749     Bool_t match = 0;
750     for (Int_t i=0;i<arr->GetEntriesFast();++i) {
751       TObject *obj = arr->At(i);
752       if (!obj)
753         continue;
754       if (fired.Contains(obj->GetName())) {
755         match = 1;
756         break;
757       }
758     }
759     delete arr;
760     if (!match) {
761       if (fGeneralHistograms) 
762         fHistEventRejection->Fill("trigger",1);
763       return kFALSE;
764     }
765   }
766
767   if (fTriggerTypeSel != kND) {
768     if (fTriggerType != fTriggerTypeSel) {
769       if (fGeneralHistograms) 
770         fHistEventRejection->Fill("trigTypeSel",1);
771       return kFALSE;
772     }
773   }
774   
775
776   if ((fMinCent != -999) && (fMaxCent != -999)) {
777     if (fCent<fMinCent || fCent>fMaxCent) {
778       if (fGeneralHistograms) 
779         fHistEventRejection->Fill("Cent",1);
780       return kFALSE;
781     }
782   }
783
784   if ((fMinVz != -999) && (fMaxVz != -999)) {
785     if (fNVertCont == 0 ) {
786       if (fGeneralHistograms) 
787         fHistEventRejection->Fill("vertex contr.",1);
788       return kFALSE;
789     }
790     Double_t vz = fVertex[2];
791     if (vz<fMinVz || vz>fMaxVz) {
792       if (fGeneralHistograms) 
793         fHistEventRejection->Fill("Vz",1);
794       return kFALSE;
795     }
796   }
797
798   if (fMinPtTrackInEmcal > 0 && fGeom) {
799     Bool_t trackInEmcalOk = kFALSE;
800     Int_t ntracks = GetNParticles(0);
801     for (Int_t i = 0; i < ntracks; i++) {
802       AliVParticle *track = GetAcceptParticleFromArray(i,0);
803       if (!track)
804         continue;
805
806       Double_t phiMin = fGeom->GetArm1PhiMin() * TMath::DegToRad();
807       Double_t phiMax = fGeom->GetArm1PhiMax() * TMath::DegToRad();
808       Int_t runNumber = InputEvent()->GetRunNumber();
809       if (runNumber>=177295 && runNumber<=197470) { //small SM masked in 2012 and 2013
810         phiMin = 1.4;   
811         phiMax = TMath::Pi();
812       }
813
814       if (track->Eta() < fGeom->GetArm1EtaMin() || track->Eta() > fGeom->GetArm1EtaMax() || track->Phi() < phiMin || track->Phi() > phiMax)
815         continue;
816       if (track->Pt() > fMinPtTrackInEmcal) {
817         trackInEmcalOk = kTRUE;
818         break;
819       }
820     }
821     if (!trackInEmcalOk) {
822       if (fGeneralHistograms) 
823         fHistEventRejection->Fill("trackInEmcal",1);
824       return kFALSE;
825     }
826   }
827
828   if (fMinNTrack > 0) {
829     Int_t nTracksAcc = 0;
830     Int_t ntracks = GetNParticles(0);
831     for (Int_t i = 0; i < ntracks; i++) {
832       AliVParticle *track = GetAcceptParticleFromArray(i,0);
833       if (!track)
834         continue;
835       if (track->Pt() > fTrackPtCut) {
836         nTracksAcc++;
837         if (nTracksAcc>=fMinNTrack)
838           break;
839       }
840     }
841     if (nTracksAcc<fMinNTrack) {
842       if (fGeneralHistograms) 
843         fHistEventRejection->Fill("minNTrack",1);
844       return kFALSE;
845     }
846   }
847
848   if (fUseAliAnaUtils) {
849     if (!fAliAnalysisUtils)
850       fAliAnalysisUtils = new AliAnalysisUtils();
851     fAliAnalysisUtils->SetMinVtxContr(2);
852     fAliAnalysisUtils->SetMaxVtxZ(10.);
853
854     if (!fAliAnalysisUtils->IsVertexSelected2013pA(InputEvent())) {
855       if (fGeneralHistograms) 
856         fHistEventRejection->Fill("VtxSel2013pA",1);
857       return kFALSE;
858     }
859
860     if (fAliAnalysisUtils->IsPileUpEvent(InputEvent())) {
861       fHistEventRejection->Fill("PileUp",1);
862       return kFALSE;
863     }
864   }
865
866   if (!(fEPV0 > fMinEventPlane && fEPV0 <= fMaxEventPlane) &&
867       !(fEPV0 + TMath::Pi() > fMinEventPlane && fEPV0 + TMath::Pi() <= fMaxEventPlane) &&
868       !(fEPV0 - TMath::Pi() > fMinEventPlane && fEPV0 - TMath::Pi() <= fMaxEventPlane)) 
869     {
870       if (fGeneralHistograms) 
871         fHistEventRejection->Fill("EvtPlane",1);
872       return kFALSE;
873     }
874
875   if (fSelectPtHardBin != -999 && fSelectPtHardBin != fPtHardBin)  {
876       if (fGeneralHistograms) 
877         fHistEventRejection->Fill("SelPtHardBin",1);
878       return kFALSE;
879     }
880
881   return kTRUE;
882 }
883
884 //________________________________________________________________________
885 TClonesArray *AliAnalysisTaskEmcal::GetArrayFromEvent(const char *name, const char *clname)
886 {
887   // Get array from event.
888
889   TClonesArray *arr = 0;
890   TString sname(name);
891   if (!sname.IsNull()) {
892     arr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(sname));
893     if (!arr) {
894       AliWarning(Form("%s: Could not retrieve array with name %s!", GetName(), name)); 
895       return 0;
896     }
897   } else {
898     return 0;
899   }
900
901   if (!clname)
902     return arr;
903
904   TString objname(arr->GetClass()->GetName());
905   TClass cls(objname);
906   if (!cls.InheritsFrom(clname)) {
907     AliWarning(Form("%s: Objects of type %s in %s are not inherited from %s!", 
908                     GetName(), cls.GetName(), name, clname)); 
909     return 0;
910   }
911   return arr;
912 }
913
914 //________________________________________________________________________
915 Bool_t AliAnalysisTaskEmcal::RetrieveEventObjects()
916 {
917   // Retrieve objects from event.
918
919   fVertex[0] = 0;
920   fVertex[1] = 0;
921   fVertex[2] = 0;
922   fNVertCont = 0;
923
924   const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
925   if (vert) {
926     vert->GetXYZ(fVertex);
927     fNVertCont = vert->GetNContributors();
928   }
929
930   fBeamType = GetBeamType();
931
932   if (fBeamType == kAA || fBeamType == kpA ) {
933     AliCentrality *aliCent = InputEvent()->GetCentrality();
934     if (aliCent) {
935       fCent = aliCent->GetCentralityPercentile(fCentEst.Data());
936       if (fNcentBins==4) {
937         if      (fCent >=  0 && fCent <   10) fCentBin = 0;
938         else if (fCent >= 10 && fCent <   30) fCentBin = 1;
939         else if (fCent >= 30 && fCent <   50) fCentBin = 2;
940         else if (fCent >= 50 && fCent <= 100) fCentBin = 3; 
941         else {
942           AliWarning(Form("%s: Negative centrality: %f. Assuming 99", GetName(), fCent));
943           fCentBin = fNcentBins-1;
944         }
945       }
946       else {
947         Double_t centWidth = (fMaxCent-fMinCent)/(Double_t)fNcentBins;
948         fCentBin = TMath::FloorNint(fCent/centWidth);
949         if (fCentBin>=fNcentBins) {
950           AliWarning(Form("%s: fCentBin too large: cent = %f fCentBin = %d. Assuming 99", GetName(),fCent,fCentBin));
951           fCentBin = fNcentBins-1;
952         }
953       }
954     } else {
955       AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
956       fCentBin = 3;
957     }
958     AliEventplane *aliEP = InputEvent()->GetEventplane();
959     if (aliEP) {
960       fEPV0  = aliEP->GetEventplane("V0" ,InputEvent());
961       fEPV0A = aliEP->GetEventplane("V0A",InputEvent());
962       fEPV0C = aliEP->GetEventplane("V0C",InputEvent());
963     } else {
964       AliWarning(Form("%s: Could not retrieve event plane information!", GetName()));
965     }
966   } else {
967     fCent = 99;
968     fCentBin = 0;
969   }
970
971   if (fIsPythia) {
972
973     if (MCEvent()) {
974       fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
975       if (!fPythiaHeader) {
976         // Check if AOD
977         AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
978
979         if (aodMCH) {
980           for (UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++) {
981             fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
982             if (fPythiaHeader) break;
983           }
984         }
985       }
986     }
987
988     if (fPythiaHeader) {
989       fPtHard = fPythiaHeader->GetPtHard();
990     
991       const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
992       const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
993       for (fPtHardBin = 0; fPtHardBin < 11; fPtHardBin++) {
994         if (fPtHard >= ptHardLo[fPtHardBin] && fPtHard < ptHardHi[fPtHardBin])
995           break;
996       }
997     
998       fNTrials = fPythiaHeader->Trials();
999     }
1000   }
1001
1002   fTriggerType = GetTriggerType();
1003
1004   return kTRUE;
1005 }
1006
1007 //________________________________________________________________________
1008 AliParticleContainer* AliAnalysisTaskEmcal::AddParticleContainer(const char *n) {
1009
1010   // Add particle container
1011   // will be called in AddTask macro
1012
1013   TString tmp = TString(n);
1014   if (tmp.IsNull()) return 0;
1015
1016   AliParticleContainer *cont = 0x0;
1017   cont = new AliParticleContainer();
1018   cont->SetArrayName(n);
1019   TString contName = cont->GetArrayName();
1020  
1021   fParticleCollArray.Add(cont);
1022
1023   return cont;
1024 }
1025
1026 //________________________________________________________________________
1027 AliClusterContainer* AliAnalysisTaskEmcal::AddClusterContainer(const char *n) {
1028
1029   // Add cluster container
1030   // will be called in AddTask macro
1031
1032   TString tmp = TString(n);
1033   if (tmp.IsNull()) return 0;
1034
1035   AliClusterContainer *cont = 0x0;
1036   cont = new AliClusterContainer();
1037   cont->SetArrayName(n);
1038
1039   fClusterCollArray.Add(cont);
1040
1041   return cont;
1042 }
1043
1044 //________________________________________________________________________
1045 AliParticleContainer* AliAnalysisTaskEmcal::GetParticleContainer(Int_t i) const {
1046   // Get i^th particle container
1047
1048   if (i<0 || i>fParticleCollArray.GetEntriesFast()) return 0;
1049   AliParticleContainer *cont = static_cast<AliParticleContainer*>(fParticleCollArray.At(i));
1050   return cont;
1051 }
1052
1053 //________________________________________________________________________
1054 AliClusterContainer* AliAnalysisTaskEmcal::GetClusterContainer(Int_t i) const {
1055   // Get i^th cluster container
1056
1057   if (i<0 || i>fClusterCollArray.GetEntriesFast()) return 0;
1058   AliClusterContainer *cont = static_cast<AliClusterContainer*>(fClusterCollArray.At(i));
1059   return cont;
1060 }
1061
1062 //________________________________________________________________________
1063 AliParticleContainer* AliAnalysisTaskEmcal::GetParticleContainer(const char *name) const {
1064   // Get particle container with name
1065
1066   AliParticleContainer *cont = static_cast<AliParticleContainer*>(fParticleCollArray.FindObject(name));
1067   return cont;
1068 }
1069
1070 //________________________________________________________________________
1071 AliClusterContainer* AliAnalysisTaskEmcal::GetClusterContainer(const char *name) const {
1072   // Get cluster container with name
1073
1074   AliClusterContainer *cont = static_cast<AliClusterContainer*>(fClusterCollArray.FindObject(name));
1075   return cont;
1076 }
1077
1078 //________________________________________________________________________
1079 TClonesArray* AliAnalysisTaskEmcal::GetParticleArray(Int_t i) const {
1080   // Get i^th TClonesArray with AliVParticle
1081
1082   AliParticleContainer *cont = GetParticleContainer(i);
1083   if (!cont) {
1084     AliError(Form("%s: Particle container %d not found",GetName(),i));
1085     return 0;
1086   }
1087   TString contName = cont->GetArrayName();
1088   return cont->GetArray();
1089 }
1090
1091 //________________________________________________________________________
1092 TClonesArray* AliAnalysisTaskEmcal::GetClusterArray(Int_t i) const {
1093   // Get i^th TClonesArray with AliVCluster
1094
1095   AliClusterContainer *cont = GetClusterContainer(i);
1096   if (!cont) {
1097     AliError(Form("%s:Cluster container %d not found",GetName(),i));
1098     return 0;
1099   }
1100   return cont->GetArray();
1101 }
1102
1103 //________________________________________________________________________
1104 AliVParticle* AliAnalysisTaskEmcal::GetAcceptParticleFromArray(Int_t p, Int_t c) const {
1105   // Get particle p if accepted from  container c
1106   // If particle not accepted return 0
1107
1108   AliParticleContainer *cont = GetParticleContainer(c);
1109   if (!cont) {
1110     AliError(Form("%s: Particle container %d not found",GetName(),c));
1111     return 0;
1112   }
1113   AliVParticle *vp = cont->GetAcceptParticle(p);
1114
1115   return vp;
1116 }
1117
1118 //________________________________________________________________________
1119 AliVCluster* AliAnalysisTaskEmcal::GetAcceptClusterFromArray(Int_t cl, Int_t c) const {
1120   // Get particle p if accepted from  container c
1121   // If particle not accepted return 0
1122
1123   AliClusterContainer *cont = GetClusterContainer(c);
1124   if (!cont) {
1125     AliError(Form("%s: Cluster container %d not found",GetName(),c));
1126     return 0;
1127   }
1128   AliVCluster *vc = cont->GetAcceptCluster(cl);
1129
1130   return vc;
1131 }
1132
1133 //________________________________________________________________________
1134 Int_t AliAnalysisTaskEmcal::GetNParticles(Int_t i) const {
1135   // Get number of entries in particle array i
1136
1137   AliParticleContainer *cont = GetParticleContainer(i);
1138   if (!cont) {
1139     AliError(Form("%s: Particle container %d not found",GetName(),i));
1140     return 0;
1141   }
1142   return cont->GetNEntries();
1143 }
1144
1145 //________________________________________________________________________
1146 Int_t AliAnalysisTaskEmcal::GetNClusters(Int_t i) const {
1147   // Get number of entries in cluster array i
1148
1149   AliClusterContainer *cont = GetClusterContainer(i);
1150   if (!cont) {
1151     AliError(Form("%s: Cluster container %d not found",GetName(),i));
1152     return 0;
1153   }
1154   return cont->GetNEntries();
1155 }
1156
1157 //________________________________________________________________________
1158 AliEmcalTriggerPatchInfo* AliAnalysisTaskEmcal::GetMainTriggerPatch() {
1159   //get main trigger match; if not known yet, look for it and cache
1160
1161   if (fMainTriggerPatch) 
1162     return fMainTriggerPatch;
1163
1164   if (!fTriggerPatchInfo) {
1165     AliError(Form("%s: fTriggerPatchInfo not available",GetName()));
1166     return 0;
1167   }
1168
1169   //number of patches in event
1170   Int_t nPatch = fTriggerPatchInfo->GetEntries();
1171
1172   //extract main trigger patch
1173   AliEmcalTriggerPatchInfo *patch;
1174   for (Int_t iPatch = 0; iPatch < nPatch; iPatch++) {
1175     
1176     patch = (AliEmcalTriggerPatchInfo*)fTriggerPatchInfo->At( iPatch );
1177     if (patch->IsMainTrigger()) {
1178       fMainTriggerPatch = patch;
1179       break;
1180     }
1181   }
1182
1183   return fMainTriggerPatch;
1184 }