3 // Emcal base analysis task.
5 // Author: S.Aiola, M. Verweij
6 #include "AliAnalysisTaskEmcal.h"
8 #include <TClonesArray.h>
18 #include "AliAODEvent.h"
19 #include "AliAnalysisManager.h"
20 #include "AliCentrality.h"
21 #include "AliEMCALGeometry.h"
22 #include "AliESDEvent.h"
23 #include "AliEmcalParticle.h"
24 #include "AliEventplane.h"
25 #include "AliInputEventHandler.h"
27 #include "AliMCParticle.h"
28 #include "AliVCluster.h"
29 #include "AliVEventHandler.h"
30 #include "AliVParticle.h"
31 #include "AliVCaloTrigger.h"
32 #include "AliGenPythiaEventHeader.h"
33 #include "AliAODMCHeader.h"
34 #include "AliMCEvent.h"
35 #include "AliAnalysisUtils.h"
36 #include "AliEmcalTriggerPatchInfo.h"
38 #include "AliParticleContainer.h"
39 #include "AliClusterContainer.h"
41 ClassImp(AliAnalysisTaskEmcal)
43 //________________________________________________________________________
44 AliAnalysisTaskEmcal::AliAnalysisTaskEmcal() :
45 AliAnalysisTaskSE("AliAnalysisTaskEmcal"),
47 fGeneralHistograms(kFALSE),
52 fCaloTriggerPatchInfoName(),
59 fUseAliAnaUtils(kFALSE),
60 fAliAnalysisUtils(0x0),
61 fOffTrigger(AliVEvent::kAny),
67 fMinPtTrackInEmcal(0),
68 fEventPlaneVsEmcal(-1),
74 fSelectPtHardBin(-999),
78 fNeedEmcalGeom(kTRUE),
100 fMainTriggerPatch(0x0),
103 fHistTrialsAfterSel(0),
104 fHistEventsAfterSel(0),
105 fHistXsectionAfterSel(0),
113 fHistEventRejection(0)
115 // Default constructor.
121 fParticleCollArray.SetOwner(kTRUE);
122 fClusterCollArray.SetOwner(kTRUE);
125 //________________________________________________________________________
126 AliAnalysisTaskEmcal::AliAnalysisTaskEmcal(const char *name, Bool_t histo) :
127 AliAnalysisTaskSE(name),
129 fGeneralHistograms(kFALSE),
130 fInitialized(kFALSE),
134 fCaloTriggerPatchInfoName(),
141 fUseAliAnaUtils(kFALSE),
142 fAliAnalysisUtils(0x0),
143 fOffTrigger(AliVEvent::kAny),
145 fTriggerTypeSel(kND),
149 fMinPtTrackInEmcal(0),
150 fEventPlaneVsEmcal(-1),
151 fMinEventPlane(-1e6),
156 fSelectPtHardBin(-999),
160 fNeedEmcalGeom(kTRUE),
167 fTriggerPatchInfo(0),
180 fParticleCollArray(),
182 fMainTriggerPatch(0x0),
185 fHistTrialsAfterSel(0),
186 fHistEventsAfterSel(0),
187 fHistXsectionAfterSel(0),
195 fHistEventRejection(0)
197 // Standard constructor.
203 fParticleCollArray.SetOwner(kTRUE);
204 fClusterCollArray.SetOwner(kTRUE);
207 DefineOutput(1, TList::Class());
211 //________________________________________________________________________
212 AliAnalysisTaskEmcal::~AliAnalysisTaskEmcal()
217 //________________________________________________________________________
218 void AliAnalysisTaskEmcal::SetClusPtCut(Double_t cut, Int_t c)
220 AliClusterContainer *cont = GetClusterContainer(c);
221 if (cont) cont->SetClusPtCut(cut);
222 else AliError(Form("%s in SetClusPtCut(...): container %d not found",GetName(),c));
225 //________________________________________________________________________
226 void AliAnalysisTaskEmcal::SetClusTimeCut(Double_t min, Double_t max, Int_t c)
228 AliClusterContainer *cont = GetClusterContainer(c);
229 if (cont) cont->SetClusTimeCut(min,max);
230 else AliError(Form("%s in SetClusTimeCut(...): container %d not found",GetName(),c));
233 //________________________________________________________________________
234 void AliAnalysisTaskEmcal::SetTrackPtCut(Double_t cut, Int_t c)
236 AliParticleContainer *cont = GetParticleContainer(c);
237 if (cont) cont->SetParticlePtCut(cut);
238 else AliError(Form("%s in SetTrackPtCut(...): container %d not found",GetName(),c));
243 //________________________________________________________________________
244 void AliAnalysisTaskEmcal::SetTrackEtaLimits(Double_t min, Double_t max, Int_t c)
246 AliParticleContainer *cont = GetParticleContainer(c);
247 if (cont) cont->SetParticleEtaLimits(min,max);
248 else AliError(Form("%s in SetTrackPtCut(...): container %d not found",GetName(),c));
251 //________________________________________________________________________
252 void AliAnalysisTaskEmcal::SetTrackPhiLimits(Double_t min, Double_t max, Int_t c)
254 AliParticleContainer *cont = GetParticleContainer(c);
255 if (cont) cont->SetParticlePhiLimits(min,max);
256 else AliError(Form("%s in SetTrackPhiLimits(...): container %d not found",GetName(),c));
259 //________________________________________________________________________
260 void AliAnalysisTaskEmcal::UserCreateOutputObjects()
262 // Create user output.
264 AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
266 AliVEventHandler *evhand = mgr->GetInputEventHandler();
268 if (evhand->InheritsFrom("AliESDInputHandler")) {
276 AliError("Event handler not found!");
280 AliError("Analysis manager not found!");
287 fOutput = new TList();
290 if (fForceBeamType == kpp)
293 if (!fGeneralHistograms)
297 fHistTrialsAfterSel = new TH1F("fHistTrialsAfterSel", "fHistTrialsAfterSel", 11, 0, 11);
298 fHistTrialsAfterSel->GetXaxis()->SetTitle("p_{T} hard bin");
299 fHistTrialsAfterSel->GetYaxis()->SetTitle("trials");
300 fOutput->Add(fHistTrialsAfterSel);
302 fHistEventsAfterSel = new TH1F("fHistEventsAfterSel", "fHistEventsAfterSel", 11, 0, 11);
303 fHistEventsAfterSel->GetXaxis()->SetTitle("p_{T} hard bin");
304 fHistEventsAfterSel->GetYaxis()->SetTitle("total events");
305 fOutput->Add(fHistEventsAfterSel);
307 fHistXsectionAfterSel = new TProfile("fHistXsectionAfterSel", "fHistXsectionAfterSel", 11, 0, 11);
308 fHistXsectionAfterSel->GetXaxis()->SetTitle("p_{T} hard bin");
309 fHistXsectionAfterSel->GetYaxis()->SetTitle("xsection");
310 fOutput->Add(fHistXsectionAfterSel);
312 fHistTrials = new TH1F("fHistTrials", "fHistTrials", 11, 0, 11);
313 fHistTrials->GetXaxis()->SetTitle("p_{T} hard bin");
314 fHistTrials->GetYaxis()->SetTitle("trials");
315 fOutput->Add(fHistTrials);
317 fHistEvents = new TH1F("fHistEvents", "fHistEvents", 11, 0, 11);
318 fHistEvents->GetXaxis()->SetTitle("p_{T} hard bin");
319 fHistEvents->GetYaxis()->SetTitle("total events");
320 fOutput->Add(fHistEvents);
322 fHistXsection = new TProfile("fHistXsection", "fHistXsection", 11, 0, 11);
323 fHistXsection->GetXaxis()->SetTitle("p_{T} hard bin");
324 fHistXsection->GetYaxis()->SetTitle("xsection");
325 fOutput->Add(fHistXsection);
327 const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
328 const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
330 for (Int_t i = 1; i < 12; i++) {
331 fHistTrialsAfterSel->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
332 fHistEventsAfterSel->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
334 fHistTrials->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
335 fHistXsection->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
336 fHistEvents->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
339 fHistPtHard = new TH1F("fHistPtHard", "fHistPtHard", fNbins*2, fMinBinPt, fMaxBinPt*4);
340 fHistPtHard->GetXaxis()->SetTitle("p_{T,hard} (GeV/c)");
341 fHistPtHard->GetYaxis()->SetTitle("counts");
342 fOutput->Add(fHistPtHard);
345 fHistCentrality = new TH1F("fHistCentrality","Event centrality distribution", 200, 0, 100);
346 fHistCentrality->GetXaxis()->SetTitle("Centrality (%)");
347 fHistCentrality->GetYaxis()->SetTitle("counts");
348 fOutput->Add(fHistCentrality);
350 fHistZVertex = new TH1F("fHistZVertex","Z vertex position", 60, -30, 30);
351 fHistZVertex->GetXaxis()->SetTitle("z");
352 fHistZVertex->GetYaxis()->SetTitle("counts");
353 fOutput->Add(fHistZVertex);
355 fHistEventPlane = new TH1F("fHistEventPlane","Event plane", 120, -TMath::Pi(), TMath::Pi());
356 fHistEventPlane->GetXaxis()->SetTitle("event plane");
357 fHistEventPlane->GetYaxis()->SetTitle("counts");
358 fOutput->Add(fHistEventPlane);
360 fHistEventRejection = new TH1I("fHistEventRejection","Reasons to reject event",20,0,20);
361 fHistEventRejection->Fill("PhysSel",0);
362 fHistEventRejection->Fill("trigger",0);
363 fHistEventRejection->Fill("trigTypeSel",0);
364 fHistEventRejection->Fill("Cent",0);
365 fHistEventRejection->Fill("vertex contr.",0);
366 fHistEventRejection->Fill("Vz",0);
367 fHistEventRejection->Fill("trackInEmcal",0);
368 fHistEventRejection->Fill("minNTrack",0);
369 fHistEventRejection->Fill("VtxSel2013pA",0);
370 fHistEventRejection->Fill("PileUp",0);
371 fHistEventRejection->Fill("EvtPlane",0);
372 fHistEventRejection->Fill("SelPtHardBin",0);
373 fHistEventRejection->GetYaxis()->SetTitle("counts");
374 fOutput->Add(fHistEventRejection);
376 PostData(1, fOutput);
379 //________________________________________________________________________
380 Bool_t AliAnalysisTaskEmcal::FillGeneralHistograms()
383 fHistEventsAfterSel->Fill(fPtHardBin, 1);
384 fHistTrialsAfterSel->Fill(fPtHardBin, fNTrials);
385 fHistXsectionAfterSel->Fill(fPtHardBin, fXsection);
386 fHistPtHard->Fill(fPtHard);
389 fHistCentrality->Fill(fCent);
390 fHistZVertex->Fill(fVertex[2]);
391 fHistEventPlane->Fill(fEPV0);
396 //________________________________________________________________________
397 void AliAnalysisTaskEmcal::UserExec(Option_t *)
399 // Main loop, called for each event.
401 fMainTriggerPatch = NULL;
409 if (!RetrieveEventObjects())
412 if (!IsEventSelected())
415 if (fGeneralHistograms && fCreateHisto) {
416 if (!FillGeneralHistograms())
424 if (!FillHistograms())
428 if (fCreateHisto && fOutput) {
429 // information for this iteration of the UserExec in the container
430 PostData(1, fOutput);
434 //________________________________________________________________________
435 Bool_t AliAnalysisTaskEmcal::AcceptCluster(AliVCluster *clus, Int_t c) const
437 // Return true if cluster is accepted.
442 AliClusterContainer *cont = GetClusterContainer(c);
444 AliError(Form("%s:Container %d not found",GetName(),c));
448 return cont->AcceptCluster(clus);
451 //________________________________________________________________________
452 Bool_t AliAnalysisTaskEmcal::AcceptTrack(AliVParticle *track, Int_t c) const
454 // Return true if track is accepted.
459 AliParticleContainer *cont = GetParticleContainer(c);
461 AliError(Form("%s:Container %d not found",GetName(),c));
465 return cont->AcceptParticle(track);
468 //________________________________________________________________________
469 Bool_t AliAnalysisTaskEmcal::PythiaInfoFromFile(const char* currFile, Float_t &fXsec, Float_t &fTrials, Int_t &pthard)
472 // Get the cross section and the trails either from pyxsec.root or from pysec_hists.root
473 // Get the pt hard bin from the file path
474 // This is to called in Notify and should provide the path to the AOD/ESD file
475 // (Partially copied from AliAnalysisHelperJetTasks)
477 TString file(currFile);
481 if (file.Contains(".zip#")) {
482 Ssiz_t pos1 = file.Index("root_archive",12,0,TString::kExact);
483 Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
484 Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
485 file.Replace(pos+1,pos2-pos1,"");
487 // not an archive take the basename....
488 file.ReplaceAll(gSystem->BaseName(file.Data()),"");
490 AliDebug(1,Form("File name: %s",file.Data()));
492 // Get the pt hard bin
493 TString strPthard(file);
495 strPthard.Remove(strPthard.Last('/'));
496 strPthard.Remove(strPthard.Last('/'));
497 if (strPthard.Contains("AOD")) strPthard.Remove(strPthard.Last('/'));
498 strPthard.Remove(0,strPthard.Last('/')+1);
499 if (strPthard.IsDec())
500 pthard = strPthard.Atoi();
502 AliWarning(Form("Could not extract file number from path %s", strPthard.Data()));
504 // 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
505 TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root"));
508 // next trial fetch the histgram file
509 fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
511 // not a severe condition but inciate that we have no information
514 // find the tlist we want to be independtent of the name so use the Tkey
515 TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
520 TList *list = dynamic_cast<TList*>(key->ReadObj());
525 fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
526 fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
529 } else { // no tree pyxsec.root
530 TTree *xtree = (TTree*)fxsec->Get("Xsection");
536 Double_t xsection = 0;
537 xtree->SetBranchAddress("xsection",&xsection);
538 xtree->SetBranchAddress("ntrials",&ntrials);
547 //________________________________________________________________________
548 Bool_t AliAnalysisTaskEmcal::UserNotify()
550 // Called when file changes.
552 if (!fIsPythia || !fGeneralHistograms || !fCreateHisto)
555 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
557 AliError(Form("%s - UserNotify: No current tree!",GetName()));
561 Float_t xsection = 0;
565 TFile *curfile = tree->GetCurrentFile();
567 AliError(Form("%s - UserNotify: No current file!",GetName()));
571 TChain *chain = dynamic_cast<TChain*>(tree);
572 if (chain) tree = chain->GetTree();
574 Int_t nevents = tree->GetEntriesFast();
576 PythiaInfoFromFile(curfile->GetName(), xsection, trials, pthardbin);
579 if ((pthardbin < 0) || (pthardbin > 10)) pthardbin = 0;
581 fHistTrials->Fill(pthardbin, trials);
582 fHistXsection->Fill(pthardbin, xsection);
583 fHistEvents->Fill(pthardbin, nevents);
588 //________________________________________________________________________
589 void AliAnalysisTaskEmcal::ExecOnce()
591 // Init the analysis.
594 AliError(Form("%s: Could not retrieve event! Returning!", GetName()));
598 if (fNeedEmcalGeom) {
599 fGeom = AliEMCALGeometry::GetInstance();
601 AliError(Form("%s: Can not create geometry", GetName()));
607 if (fEventPlaneVsEmcal >= 0) {
609 Double_t ep = (fGeom->GetArm1PhiMax() + fGeom->GetArm1PhiMin()) / 2 * TMath::DegToRad() + fEventPlaneVsEmcal - TMath::Pi();
610 fMinEventPlane = ep - TMath::Pi() / 4;
611 fMaxEventPlane = ep + TMath::Pi() / 4;
614 AliWarning("Could not set event plane limits because EMCal geometry was not loaded!");
618 //Load all requested track branches - each container knows name already
619 for (Int_t i =0; i<fParticleCollArray.GetEntriesFast(); i++) {
620 AliParticleContainer *cont = static_cast<AliParticleContainer*>(fParticleCollArray.At(i));
621 cont->SetArray(InputEvent());
624 if (fParticleCollArray.GetEntriesFast()>0) {
625 fTracks = GetParticleArray(0);
627 AliError(Form("%s: Could not retrieve first track branch!", GetName()));
632 //Load all requested cluster branches - each container knows name already
633 for (Int_t i =0; i<fClusterCollArray.GetEntriesFast(); i++) {
634 AliClusterContainer *cont = static_cast<AliClusterContainer*>(fClusterCollArray.At(i));
635 cont->SetArray(InputEvent());
638 if (fClusterCollArray.GetEntriesFast()>0) {
639 fCaloClusters = GetClusterArray(0);
640 if (!fCaloClusters) {
641 AliError(Form("%s: Could not retrieve first cluster branch!", GetName()));
646 if (!fCaloCellsName.IsNull() && !fCaloCells) {
647 fCaloCells = dynamic_cast<AliVCaloCells*>(InputEvent()->FindListObject(fCaloCellsName));
649 AliError(Form("%s: Could not retrieve cells %s!", GetName(), fCaloCellsName.Data()));
654 if (!fCaloTriggersName.IsNull() && !fCaloTriggers) {
655 fCaloTriggers = dynamic_cast<AliVCaloTrigger*>(InputEvent()->FindListObject(fCaloTriggersName));
656 if (!fCaloTriggers) {
657 AliError(Form("%s: Could not retrieve calo triggers %s!", GetName(), fCaloTriggersName.Data()));
662 if (!fCaloTriggerPatchInfoName.IsNull() && !fTriggerPatchInfo) {
663 fTriggerPatchInfo = GetArrayFromEvent(fCaloTriggerPatchInfoName.Data(),"AliEmcalTriggerPatchInfo");
664 if (!fTriggerPatchInfo) {
665 AliError(Form("%s: Could not retrieve calo trigger patch info %s!", GetName(), fCaloTriggerPatchInfoName.Data()));
671 fInitialized = kTRUE;
674 //_____________________________________________________
675 AliAnalysisTaskEmcal::BeamType AliAnalysisTaskEmcal::GetBeamType()
677 // Get beam type : pp-AA-pA
678 // ESDs have it directly, AODs get it from hardcoded run number ranges
680 if (fForceBeamType != kNA)
681 return fForceBeamType;
683 AliESDEvent *esd = dynamic_cast<AliESDEvent*>(InputEvent());
685 const AliESDRun *run = esd->GetESDRun();
686 TString beamType = run->GetBeamType();
687 if (beamType == "p-p")
689 else if (beamType == "A-A")
691 else if (beamType == "p-A")
696 Int_t runNumber = InputEvent()->GetRunNumber();
697 if ((runNumber >= 136851 && runNumber <= 139517) || // LHC10h
698 (runNumber >= 166529 && runNumber <= 170593)) { // LHC11h
700 } else if ((runNumber>=188365 && runNumber <= 188366) || // LHC12g
701 (runNumber >= 195344 && runNumber <= 196608)) { // LHC13b-f
709 //________________________________________________________________________
710 ULong_t AliAnalysisTaskEmcal::GetTriggerList(){
711 if (!fTriggerPatchInfo)
714 //number of patches in event
715 Int_t nPatch = fTriggerPatchInfo->GetEntries();
717 //loop over patches to define trigger type of event
722 AliEmcalTriggerPatchInfo *patch;
723 for (Int_t iPatch = 0; iPatch < nPatch; iPatch++) {
724 patch = (AliEmcalTriggerPatchInfo*)fTriggerPatchInfo->At( iPatch );
725 if (patch->IsGammaHigh()) nG1++;
726 if (patch->IsGammaLow()) nG2++;
727 if (patch->IsJetHigh()) nJ1++;
728 if (patch->IsJetLow()) nJ2++;
731 AliDebug(2, "Patch summary: ");
732 AliDebug(2, Form("Number of patches: %d", nPatch));
733 AliDebug(2, Form("Jet: low[%d], high[%d]" ,nJ2, nJ1));
734 AliDebug(2, Form("Gamma: low[%d], high[%d]" ,nG2, nG1));
738 SETBIT(triggers, kG1);
740 SETBIT(triggers, kG2);
742 SETBIT(triggers, kJ1);
744 SETBIT(triggers, kJ2);
748 //________________________________________________________________________
749 Bool_t AliAnalysisTaskEmcal::HasTriggerType(TriggerType t){
750 // Check if event has a given trigger type
752 return fTriggers == 0;
754 return TESTBIT(fTriggers, int(t));
757 //________________________________________________________________________
758 Bool_t AliAnalysisTaskEmcal::IsEventSelected()
760 // Check if event is selected
762 if (fOffTrigger != AliVEvent::kAny) {
764 const AliESDEvent *eev = dynamic_cast<const AliESDEvent*>(InputEvent());
766 res = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
768 const AliAODEvent *aev = dynamic_cast<const AliAODEvent*>(InputEvent());
770 res = aev->GetHeader()->GetOfflineTrigger();
773 if ((res & fOffTrigger) == 0) {
774 if (fGeneralHistograms)
775 fHistEventRejection->Fill("PhysSel",1);
780 if (!fTrigClass.IsNull()) {
782 const AliESDEvent *eev = dynamic_cast<const AliESDEvent*>(InputEvent());
784 fired = eev->GetFiredTriggerClasses();
786 const AliAODEvent *aev = dynamic_cast<const AliAODEvent*>(InputEvent());
788 fired = aev->GetFiredTriggerClasses();
791 if (!fired.Contains("-B-")) {
792 if (fGeneralHistograms)
793 fHistEventRejection->Fill("trigger",1);
796 TObjArray *arr = fTrigClass.Tokenize("|");
798 if (fGeneralHistograms)
799 fHistEventRejection->Fill("trigger",1);
803 for (Int_t i=0;i<arr->GetEntriesFast();++i) {
804 TObject *obj = arr->At(i);
807 if (fired.Contains(obj->GetName())) {
814 if (fGeneralHistograms)
815 fHistEventRejection->Fill("trigger",1);
820 if (fTriggerTypeSel != kND) {
821 if (!HasTriggerType(fTriggerTypeSel)) {
822 if (fGeneralHistograms)
823 fHistEventRejection->Fill("trigTypeSel",1);
828 if ((fMinCent != -999) && (fMaxCent != -999)) {
829 if (fCent<fMinCent || fCent>fMaxCent) {
830 if (fGeneralHistograms)
831 fHistEventRejection->Fill("Cent",1);
836 if (fUseAliAnaUtils) {
837 if (!fAliAnalysisUtils)
838 fAliAnalysisUtils = new AliAnalysisUtils();
839 fAliAnalysisUtils->SetMinVtxContr(2);
840 fAliAnalysisUtils->SetMaxVtxZ(999);
841 if(fMinVz<-10.) fMinVz = -10.;
842 if(fMinVz>10.) fMaxVz = 10.;
844 if (!fAliAnalysisUtils->IsVertexSelected2013pA(InputEvent())) {
845 if (fGeneralHistograms)
846 fHistEventRejection->Fill("VtxSel2013pA",1);
850 if (fAliAnalysisUtils->IsPileUpEvent(InputEvent())) {
851 fHistEventRejection->Fill("PileUp",1);
856 if ((fMinVz != -999) && (fMaxVz != -999)) {
857 if (fNVertCont == 0 ) {
858 if (fGeneralHistograms)
859 fHistEventRejection->Fill("vertex contr.",1);
862 Double_t vz = fVertex[2];
863 if (vz<fMinVz || vz>fMaxVz) {
864 if (fGeneralHistograms)
865 fHistEventRejection->Fill("Vz",1);
870 if (fMinPtTrackInEmcal > 0 && fGeom) {
871 Bool_t trackInEmcalOk = kFALSE;
872 Int_t ntracks = GetNParticles(0);
873 for (Int_t i = 0; i < ntracks; i++) {
874 AliVParticle *track = GetAcceptParticleFromArray(i,0);
878 Double_t phiMin = fGeom->GetArm1PhiMin() * TMath::DegToRad();
879 Double_t phiMax = fGeom->GetArm1PhiMax() * TMath::DegToRad();
880 Int_t runNumber = InputEvent()->GetRunNumber();
881 if (runNumber>=177295 && runNumber<=197470) { //small SM masked in 2012 and 2013
883 phiMax = TMath::Pi();
886 if (track->Eta() < fGeom->GetArm1EtaMin() || track->Eta() > fGeom->GetArm1EtaMax() || track->Phi() < phiMin || track->Phi() > phiMax)
888 if (track->Pt() > fMinPtTrackInEmcal) {
889 trackInEmcalOk = kTRUE;
893 if (!trackInEmcalOk) {
894 if (fGeneralHistograms)
895 fHistEventRejection->Fill("trackInEmcal",1);
900 if (fMinNTrack > 0) {
901 Int_t nTracksAcc = 0;
902 Int_t ntracks = GetNParticles(0);
903 for (Int_t i = 0; i < ntracks; i++) {
904 AliVParticle *track = GetAcceptParticleFromArray(i,0);
907 if (track->Pt() > fTrackPtCut) {
909 if (nTracksAcc>=fMinNTrack)
913 if (nTracksAcc<fMinNTrack) {
914 if (fGeneralHistograms)
915 fHistEventRejection->Fill("minNTrack",1);
920 if (!(fEPV0 > fMinEventPlane && fEPV0 <= fMaxEventPlane) &&
921 !(fEPV0 + TMath::Pi() > fMinEventPlane && fEPV0 + TMath::Pi() <= fMaxEventPlane) &&
922 !(fEPV0 - TMath::Pi() > fMinEventPlane && fEPV0 - TMath::Pi() <= fMaxEventPlane))
924 if (fGeneralHistograms)
925 fHistEventRejection->Fill("EvtPlane",1);
929 if (fSelectPtHardBin != -999 && fSelectPtHardBin != fPtHardBin) {
930 if (fGeneralHistograms)
931 fHistEventRejection->Fill("SelPtHardBin",1);
938 //________________________________________________________________________
939 TClonesArray *AliAnalysisTaskEmcal::GetArrayFromEvent(const char *name, const char *clname)
941 // Get array from event.
943 TClonesArray *arr = 0;
945 if (!sname.IsNull()) {
946 arr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(sname));
948 AliWarning(Form("%s: Could not retrieve array with name %s!", GetName(), name));
958 TString objname(arr->GetClass()->GetName());
960 if (!cls.InheritsFrom(clname)) {
961 AliWarning(Form("%s: Objects of type %s in %s are not inherited from %s!",
962 GetName(), cls.GetName(), name, clname));
968 //________________________________________________________________________
969 Bool_t AliAnalysisTaskEmcal::RetrieveEventObjects()
971 // Retrieve objects from event.
978 const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
980 vert->GetXYZ(fVertex);
981 fNVertCont = vert->GetNContributors();
984 fBeamType = GetBeamType();
986 if (fBeamType == kAA || fBeamType == kpA ) {
987 AliCentrality *aliCent = InputEvent()->GetCentrality();
989 fCent = aliCent->GetCentralityPercentile(fCentEst.Data());
991 if (fCent >= 0 && fCent < 10) fCentBin = 0;
992 else if (fCent >= 10 && fCent < 30) fCentBin = 1;
993 else if (fCent >= 30 && fCent < 50) fCentBin = 2;
994 else if (fCent >= 50 && fCent <= 100) fCentBin = 3;
996 AliWarning(Form("%s: Negative centrality: %f. Assuming 99", GetName(), fCent));
997 fCentBin = fNcentBins-1;
1000 Double_t centWidth = (fMaxCent-fMinCent)/(Double_t)fNcentBins;
1002 fCentBin = TMath::FloorNint(fCent/centWidth);
1005 if (fCentBin>=fNcentBins) {
1006 AliWarning(Form("%s: fCentBin too large: cent = %f fCentBin = %d. Assuming 99", GetName(),fCent,fCentBin));
1007 fCentBin = fNcentBins-1;
1011 AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
1012 fCentBin = fNcentBins-1;
1014 AliEventplane *aliEP = InputEvent()->GetEventplane();
1016 fEPV0 = aliEP->GetEventplane("V0" ,InputEvent());
1017 fEPV0A = aliEP->GetEventplane("V0A",InputEvent());
1018 fEPV0C = aliEP->GetEventplane("V0C",InputEvent());
1020 AliWarning(Form("%s: Could not retrieve event plane information!", GetName()));
1030 fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
1031 if (!fPythiaHeader) {
1033 AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
1036 for (UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++) {
1037 fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
1038 if (fPythiaHeader) break;
1044 if (fPythiaHeader) {
1045 fPtHard = fPythiaHeader->GetPtHard();
1047 const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
1048 const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
1049 for (fPtHardBin = 0; fPtHardBin < 11; fPtHardBin++) {
1050 if (fPtHard >= ptHardLo[fPtHardBin] && fPtHard < ptHardHi[fPtHardBin])
1054 fXsection = fPythiaHeader->GetXsection();
1055 fNTrials = fPythiaHeader->Trials();
1059 fTriggers = GetTriggerList();
1064 //________________________________________________________________________
1065 AliParticleContainer* AliAnalysisTaskEmcal::AddParticleContainer(const char *n)
1067 // Add particle container
1068 // will be called in AddTask macro
1070 TString tmp = TString(n);
1071 if (tmp.IsNull()) return 0;
1073 AliParticleContainer *cont = 0x0;
1074 cont = new AliParticleContainer();
1075 cont->SetArrayName(n);
1076 TString contName = cont->GetArrayName();
1078 fParticleCollArray.Add(cont);
1083 //________________________________________________________________________
1084 AliClusterContainer* AliAnalysisTaskEmcal::AddClusterContainer(const char *n)
1086 // Add cluster container
1087 // will be called in AddTask macro
1089 TString tmp = TString(n);
1090 if (tmp.IsNull()) return 0;
1092 AliClusterContainer *cont = 0x0;
1093 cont = new AliClusterContainer();
1094 cont->SetArrayName(n);
1096 fClusterCollArray.Add(cont);
1101 //________________________________________________________________________
1102 AliParticleContainer* AliAnalysisTaskEmcal::GetParticleContainer(Int_t i) const
1104 // Get i^th particle container
1106 if (i<0 || i>fParticleCollArray.GetEntriesFast()) return 0;
1107 AliParticleContainer *cont = static_cast<AliParticleContainer*>(fParticleCollArray.At(i));
1111 //________________________________________________________________________
1112 AliClusterContainer* AliAnalysisTaskEmcal::GetClusterContainer(Int_t i) const
1114 // Get i^th cluster container
1116 if (i<0 || i>fClusterCollArray.GetEntriesFast()) return 0;
1117 AliClusterContainer *cont = static_cast<AliClusterContainer*>(fClusterCollArray.At(i));
1121 //________________________________________________________________________
1122 AliParticleContainer* AliAnalysisTaskEmcal::GetParticleContainer(const char *name) const
1124 // Get particle container with name
1126 AliParticleContainer *cont = static_cast<AliParticleContainer*>(fParticleCollArray.FindObject(name));
1130 //________________________________________________________________________
1131 AliClusterContainer* AliAnalysisTaskEmcal::GetClusterContainer(const char *name) const
1133 // Get cluster container with name
1135 AliClusterContainer *cont = static_cast<AliClusterContainer*>(fClusterCollArray.FindObject(name));
1139 //________________________________________________________________________
1140 TClonesArray* AliAnalysisTaskEmcal::GetParticleArray(Int_t i) const
1142 // Get i^th TClonesArray with AliVParticle
1144 AliParticleContainer *cont = GetParticleContainer(i);
1146 AliError(Form("%s: Particle container %d not found",GetName(),i));
1149 TString contName = cont->GetArrayName();
1150 return cont->GetArray();
1153 //________________________________________________________________________
1154 TClonesArray* AliAnalysisTaskEmcal::GetClusterArray(Int_t i) const
1156 // Get i^th TClonesArray with AliVCluster
1158 AliClusterContainer *cont = GetClusterContainer(i);
1160 AliError(Form("%s:Cluster container %d not found",GetName(),i));
1163 return cont->GetArray();
1166 //________________________________________________________________________
1167 AliVParticle* AliAnalysisTaskEmcal::GetAcceptParticleFromArray(Int_t p, Int_t c) const
1169 // Get particle p if accepted from container c
1170 // If particle not accepted return 0
1172 AliParticleContainer *cont = GetParticleContainer(c);
1174 AliError(Form("%s: Particle container %d not found",GetName(),c));
1177 AliVParticle *vp = cont->GetAcceptParticle(p);
1182 //________________________________________________________________________
1183 AliVCluster* AliAnalysisTaskEmcal::GetAcceptClusterFromArray(Int_t cl, Int_t c) const
1185 // Get particle p if accepted from container c
1186 // If particle not accepted return 0
1188 AliClusterContainer *cont = GetClusterContainer(c);
1190 AliError(Form("%s: Cluster container %d not found",GetName(),c));
1193 AliVCluster *vc = cont->GetAcceptCluster(cl);
1198 //________________________________________________________________________
1199 Int_t AliAnalysisTaskEmcal::GetNParticles(Int_t i) const
1201 // Get number of entries in particle array i
1203 AliParticleContainer *cont = GetParticleContainer(i);
1205 AliError(Form("%s: Particle container %d not found",GetName(),i));
1208 return cont->GetNEntries();
1211 //________________________________________________________________________
1212 Int_t AliAnalysisTaskEmcal::GetNClusters(Int_t i) const
1214 // Get number of entries in cluster array i
1216 AliClusterContainer *cont = GetClusterContainer(i);
1218 AliError(Form("%s: Cluster container %d not found",GetName(),i));
1221 return cont->GetNEntries();
1224 //________________________________________________________________________
1225 AliEmcalTriggerPatchInfo* AliAnalysisTaskEmcal::GetMainTriggerPatch()
1227 //get main trigger match; if not known yet, look for it and cache
1229 if (fMainTriggerPatch)
1230 return fMainTriggerPatch;
1232 if (!fTriggerPatchInfo) {
1233 AliError(Form("%s: fTriggerPatchInfo not available",GetName()));
1237 //number of patches in event
1238 Int_t nPatch = fTriggerPatchInfo->GetEntries();
1240 //extract main trigger patch
1241 AliEmcalTriggerPatchInfo *patch;
1242 for (Int_t iPatch = 0; iPatch < nPatch; iPatch++) {
1244 patch = (AliEmcalTriggerPatchInfo*)fTriggerPatchInfo->At( iPatch );
1245 if (patch->IsMainTrigger()) {
1246 fMainTriggerPatch = patch;
1251 return fMainTriggerPatch;
1254 //________________________________________________________________________
1255 void AliAnalysisTaskEmcal::AddObjectToEvent(TObject *obj)
1257 // Add object to event
1259 if (!(InputEvent()->FindListObject(obj->GetName()))) {
1260 InputEvent()->AddObject(obj);
1262 AliFatal(Form("%s: Container with name %s already present. Aborting", GetName(), obj->GetName()));
1266 //________________________________________________________________________
1267 void AliAnalysisTaskEmcal::GenerateFixedBinArray(Int_t n, Double_t min, Double_t max, Double_t* array) const
1269 Double_t binWidth = (max-min)/n;
1271 for (Int_t i = 1; i <= n; i++) {
1272 array[i] = array[i-1]+binWidth;
1276 //________________________________________________________________________
1277 Double_t* AliAnalysisTaskEmcal::GenerateFixedBinArray(Int_t n, Double_t min, Double_t max) const
1279 Double_t *array = new Double_t[n+1];
1281 GenerateFixedBinArray(n, min, max, array);