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