2 // Emcal base analysis task.
4 // 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),
50 fMainTriggerPatchSet(kFALSE),
53 fCaloTriggerPatchInfoName(),
60 fUseAliAnaUtils(kFALSE),
61 fRejectPileup(kFALSE),
62 fTklVsClusSPDCut(kFALSE),
63 fAliAnalysisUtils(0x0),
64 fOffTrigger(AliVEvent::kAny),
70 fMinPtTrackInEmcal(0),
71 fEventPlaneVsEmcal(-1),
77 fSelectPtHardBin(-999),
81 fNeedEmcalGeom(kTRUE),
101 fParticleCollArray(),
106 fHistTrialsAfterSel(0),
107 fHistEventsAfterSel(0),
108 fHistXsectionAfterSel(0),
116 fHistEventRejection(0)
118 // Default constructor.
124 fParticleCollArray.SetOwner(kTRUE);
125 fClusterCollArray.SetOwner(kTRUE);
126 memset(fMainTriggerPatch, 0, sizeof(AliEmcalTriggerPatchInfo *) * kNType);
129 //________________________________________________________________________
130 AliAnalysisTaskEmcal::AliAnalysisTaskEmcal(const char *name, Bool_t histo) :
131 AliAnalysisTaskSE(name),
133 fGeneralHistograms(kFALSE),
134 fInitialized(kFALSE),
136 fMainTriggerPatchSet(kFALSE),
139 fCaloTriggerPatchInfoName(),
146 fUseAliAnaUtils(kFALSE),
147 fRejectPileup(kFALSE),
148 fTklVsClusSPDCut(kFALSE),
149 fAliAnalysisUtils(0x0),
150 fOffTrigger(AliVEvent::kAny),
152 fTriggerTypeSel(kND),
156 fMinPtTrackInEmcal(0),
157 fEventPlaneVsEmcal(-1),
158 fMinEventPlane(-1e6),
163 fSelectPtHardBin(-999),
167 fNeedEmcalGeom(kTRUE),
174 fTriggerPatchInfo(0),
187 fParticleCollArray(),
192 fHistTrialsAfterSel(0),
193 fHistEventsAfterSel(0),
194 fHistXsectionAfterSel(0),
202 fHistEventRejection(0)
204 // Standard constructor.
210 fParticleCollArray.SetOwner(kTRUE);
211 fClusterCollArray.SetOwner(kTRUE);
212 memset(fMainTriggerPatch, 0, sizeof(AliEmcalTriggerPatchInfo *) * kNType);
215 DefineOutput(1, TList::Class());
219 //________________________________________________________________________
220 AliAnalysisTaskEmcal::~AliAnalysisTaskEmcal()
225 //________________________________________________________________________
226 void AliAnalysisTaskEmcal::SetClusPtCut(Double_t cut, Int_t c)
228 AliClusterContainer *cont = GetClusterContainer(c);
229 if (cont) cont->SetClusPtCut(cut);
230 else AliError(Form("%s in SetClusPtCut(...): container %d not found",GetName(),c));
233 //________________________________________________________________________
234 void AliAnalysisTaskEmcal::SetClusTimeCut(Double_t min, Double_t max, Int_t c)
236 AliClusterContainer *cont = GetClusterContainer(c);
237 if (cont) cont->SetClusTimeCut(min,max);
238 else AliError(Form("%s in SetClusTimeCut(...): container %d not found",GetName(),c));
241 //________________________________________________________________________
242 void AliAnalysisTaskEmcal::SetTrackPtCut(Double_t cut, Int_t c)
244 AliParticleContainer *cont = GetParticleContainer(c);
245 if (cont) cont->SetParticlePtCut(cut);
246 else AliError(Form("%s in SetTrackPtCut(...): container %d not found",GetName(),c));
251 //________________________________________________________________________
252 void AliAnalysisTaskEmcal::SetTrackEtaLimits(Double_t min, Double_t max, Int_t c)
254 AliParticleContainer *cont = GetParticleContainer(c);
255 if (cont) cont->SetParticleEtaLimits(min,max);
256 else AliError(Form("%s in SetTrackPtCut(...): container %d not found",GetName(),c));
259 //________________________________________________________________________
260 void AliAnalysisTaskEmcal::SetTrackPhiLimits(Double_t min, Double_t max, Int_t c)
262 AliParticleContainer *cont = GetParticleContainer(c);
263 if (cont) cont->SetParticlePhiLimits(min,max);
264 else AliError(Form("%s in SetTrackPhiLimits(...): container %d not found",GetName(),c));
267 //________________________________________________________________________
268 void AliAnalysisTaskEmcal::UserCreateOutputObjects()
270 // Create user output.
272 AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
274 AliVEventHandler *evhand = mgr->GetInputEventHandler();
276 if (evhand->InheritsFrom("AliESDInputHandler")) {
284 AliError("Event handler not found!");
288 AliError("Analysis manager not found!");
295 fOutput = new TList();
298 if (fForceBeamType == kpp)
301 if (!fGeneralHistograms)
305 fHistTrialsAfterSel = new TH1F("fHistTrialsAfterSel", "fHistTrialsAfterSel", 11, 0, 11);
306 fHistTrialsAfterSel->GetXaxis()->SetTitle("p_{T} hard bin");
307 fHistTrialsAfterSel->GetYaxis()->SetTitle("trials");
308 fOutput->Add(fHistTrialsAfterSel);
310 fHistEventsAfterSel = new TH1F("fHistEventsAfterSel", "fHistEventsAfterSel", 11, 0, 11);
311 fHistEventsAfterSel->GetXaxis()->SetTitle("p_{T} hard bin");
312 fHistEventsAfterSel->GetYaxis()->SetTitle("total events");
313 fOutput->Add(fHistEventsAfterSel);
315 fHistXsectionAfterSel = new TProfile("fHistXsectionAfterSel", "fHistXsectionAfterSel", 11, 0, 11);
316 fHistXsectionAfterSel->GetXaxis()->SetTitle("p_{T} hard bin");
317 fHistXsectionAfterSel->GetYaxis()->SetTitle("xsection");
318 fOutput->Add(fHistXsectionAfterSel);
320 fHistTrials = new TH1F("fHistTrials", "fHistTrials", 11, 0, 11);
321 fHistTrials->GetXaxis()->SetTitle("p_{T} hard bin");
322 fHistTrials->GetYaxis()->SetTitle("trials");
323 fOutput->Add(fHistTrials);
325 fHistEvents = new TH1F("fHistEvents", "fHistEvents", 11, 0, 11);
326 fHistEvents->GetXaxis()->SetTitle("p_{T} hard bin");
327 fHistEvents->GetYaxis()->SetTitle("total events");
328 fOutput->Add(fHistEvents);
330 fHistXsection = new TProfile("fHistXsection", "fHistXsection", 11, 0, 11);
331 fHistXsection->GetXaxis()->SetTitle("p_{T} hard bin");
332 fHistXsection->GetYaxis()->SetTitle("xsection");
333 fOutput->Add(fHistXsection);
335 const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
336 const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
338 for (Int_t i = 1; i < 12; i++) {
339 fHistTrialsAfterSel->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
340 fHistEventsAfterSel->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
342 fHistTrials->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
343 fHistXsection->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
344 fHistEvents->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
347 fHistPtHard = new TH1F("fHistPtHard", "fHistPtHard", fNbins*2, fMinBinPt, fMaxBinPt*4);
348 fHistPtHard->GetXaxis()->SetTitle("p_{T,hard} (GeV/c)");
349 fHistPtHard->GetYaxis()->SetTitle("counts");
350 fOutput->Add(fHistPtHard);
353 fHistCentrality = new TH1F("fHistCentrality","Event centrality distribution", 200, 0, 100);
354 fHistCentrality->GetXaxis()->SetTitle("Centrality (%)");
355 fHistCentrality->GetYaxis()->SetTitle("counts");
356 fOutput->Add(fHistCentrality);
358 fHistZVertex = new TH1F("fHistZVertex","Z vertex position", 60, -30, 30);
359 fHistZVertex->GetXaxis()->SetTitle("z");
360 fHistZVertex->GetYaxis()->SetTitle("counts");
361 fOutput->Add(fHistZVertex);
363 fHistEventPlane = new TH1F("fHistEventPlane","Event plane", 120, -TMath::Pi(), TMath::Pi());
364 fHistEventPlane->GetXaxis()->SetTitle("event plane");
365 fHistEventPlane->GetYaxis()->SetTitle("counts");
366 fOutput->Add(fHistEventPlane);
368 fHistEventRejection = new TH1F("fHistEventRejection","Reasons to reject event",20,0,20);
369 fHistEventRejection->GetXaxis()->SetBinLabel(1,"PhysSel");
370 fHistEventRejection->GetXaxis()->SetBinLabel(2,"trigger");
371 fHistEventRejection->GetXaxis()->SetBinLabel(3,"trigTypeSel");
372 fHistEventRejection->GetXaxis()->SetBinLabel(4,"Cent");
373 fHistEventRejection->GetXaxis()->SetBinLabel(5,"vertex contr.");
374 fHistEventRejection->GetXaxis()->SetBinLabel(6,"Vz");
375 fHistEventRejection->GetXaxis()->SetBinLabel(7,"trackInEmcal");
376 fHistEventRejection->GetXaxis()->SetBinLabel(8,"minNTrack");
377 fHistEventRejection->GetXaxis()->SetBinLabel(9,"VtxSel2013pA");
378 fHistEventRejection->GetXaxis()->SetBinLabel(10,"PileUp");
379 fHistEventRejection->GetXaxis()->SetBinLabel(11,"EvtPlane");
380 fHistEventRejection->GetXaxis()->SetBinLabel(12,"SelPtHardBin");
381 fHistEventRejection->GetXaxis()->SetBinLabel(13,"Bkg evt");
382 fHistEventRejection->GetYaxis()->SetTitle("counts");
383 fOutput->Add(fHistEventRejection);
385 fHistEventCount = new TH1F("fHistEventCount","fHistEventCount",2,0,2);
386 fHistEventCount->GetXaxis()->SetBinLabel(1,"Accepted");
387 fHistEventCount->GetXaxis()->SetBinLabel(2,"Rejected");
388 fHistEventCount->GetYaxis()->SetTitle("counts");
389 fOutput->Add(fHistEventCount);
391 PostData(1, fOutput);
394 //________________________________________________________________________
395 Bool_t AliAnalysisTaskEmcal::FillGeneralHistograms()
398 fHistEventsAfterSel->Fill(fPtHardBin, 1);
399 fHistTrialsAfterSel->Fill(fPtHardBin, fNTrials);
400 fHistXsectionAfterSel->Fill(fPtHardBin, fXsection);
401 fHistPtHard->Fill(fPtHard);
404 fHistCentrality->Fill(fCent);
405 fHistZVertex->Fill(fVertex[2]);
406 fHistEventPlane->Fill(fEPV0);
411 //________________________________________________________________________
412 void AliAnalysisTaskEmcal::UserExec(Option_t *)
414 // Main loop, called for each event.
416 memset(fMainTriggerPatch, 0, sizeof(AliEmcalTriggerPatchInfo *) * kNType);
424 if (!RetrieveEventObjects())
427 if (IsEventSelected()) {
428 if (fGeneralHistograms) fHistEventCount->Fill("Accepted",1);
431 if (fGeneralHistograms) fHistEventCount->Fill("Rejected",1);
435 if (fGeneralHistograms && fCreateHisto) {
436 if (!FillGeneralHistograms())
444 if (!FillHistograms())
448 if (fCreateHisto && fOutput) {
449 // information for this iteration of the UserExec in the container
450 PostData(1, fOutput);
454 //________________________________________________________________________
455 Bool_t AliAnalysisTaskEmcal::AcceptCluster(AliVCluster *clus, Int_t c) const
457 // Return true if cluster is accepted.
462 AliClusterContainer *cont = GetClusterContainer(c);
464 AliError(Form("%s:Container %d not found",GetName(),c));
468 return cont->AcceptCluster(clus);
471 //________________________________________________________________________
472 Bool_t AliAnalysisTaskEmcal::AcceptTrack(AliVParticle *track, Int_t c) const
474 // Return true if track is accepted.
479 AliParticleContainer *cont = GetParticleContainer(c);
481 AliError(Form("%s:Container %d not found",GetName(),c));
485 return cont->AcceptParticle(track);
488 //________________________________________________________________________
489 Bool_t AliAnalysisTaskEmcal::PythiaInfoFromFile(const char* currFile, Float_t &fXsec, Float_t &fTrials, Int_t &pthard)
492 // Get the cross section and the trails either from pyxsec.root or from pysec_hists.root
493 // Get the pt hard bin from the file path
494 // This is to called in Notify and should provide the path to the AOD/ESD file
495 // (Partially copied from AliAnalysisHelperJetTasks)
497 TString file(currFile);
501 if (file.Contains(".zip#")) {
502 Ssiz_t pos1 = file.Index("root_archive",12,0,TString::kExact);
503 Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
504 Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
505 file.Replace(pos+1,pos2-pos1,"");
507 // not an archive take the basename....
508 file.ReplaceAll(gSystem->BaseName(file.Data()),"");
510 AliDebug(1,Form("File name: %s",file.Data()));
512 // Get the pt hard bin
513 TString strPthard(file);
515 strPthard.Remove(strPthard.Last('/'));
516 strPthard.Remove(strPthard.Last('/'));
517 if (strPthard.Contains("AOD")) strPthard.Remove(strPthard.Last('/'));
518 strPthard.Remove(0,strPthard.Last('/')+1);
519 if (strPthard.IsDec())
520 pthard = strPthard.Atoi();
522 AliWarning(Form("Could not extract file number from path %s", strPthard.Data()));
524 // 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
525 TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root"));
528 // next trial fetch the histgram file
529 fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
531 // not a severe condition but inciate that we have no information
534 // find the tlist we want to be independtent of the name so use the Tkey
535 TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
540 TList *list = dynamic_cast<TList*>(key->ReadObj());
545 fXsec = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
546 fTrials = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
549 } else { // no tree pyxsec.root
550 TTree *xtree = (TTree*)fxsec->Get("Xsection");
556 Double_t xsection = 0;
557 xtree->SetBranchAddress("xsection",&xsection);
558 xtree->SetBranchAddress("ntrials",&ntrials);
567 //________________________________________________________________________
568 Bool_t AliAnalysisTaskEmcal::UserNotify()
570 // Called when file changes.
572 if (!fIsPythia || !fGeneralHistograms || !fCreateHisto)
575 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
577 AliError(Form("%s - UserNotify: No current tree!",GetName()));
581 Float_t xsection = 0;
585 TFile *curfile = tree->GetCurrentFile();
587 AliError(Form("%s - UserNotify: No current file!",GetName()));
591 TChain *chain = dynamic_cast<TChain*>(tree);
592 if (chain) tree = chain->GetTree();
594 Int_t nevents = tree->GetEntriesFast();
596 PythiaInfoFromFile(curfile->GetName(), xsection, trials, pthardbin);
599 if ((pthardbin < 0) || (pthardbin > 10)) pthardbin = 0;
601 fHistTrials->Fill(pthardbin, trials);
602 fHistXsection->Fill(pthardbin, xsection);
603 fHistEvents->Fill(pthardbin, nevents);
608 //________________________________________________________________________
609 void AliAnalysisTaskEmcal::ExecOnce()
611 // Init the analysis.
614 AliError(Form("%s: Could not retrieve event! Returning!", GetName()));
618 if (fNeedEmcalGeom) {
619 fGeom = AliEMCALGeometry::GetInstance();
621 AliError(Form("%s: Can not create geometry", GetName()));
627 if (fEventPlaneVsEmcal >= 0) {
629 Double_t ep = (fGeom->GetArm1PhiMax() + fGeom->GetArm1PhiMin()) / 2 * TMath::DegToRad() + fEventPlaneVsEmcal - TMath::Pi();
630 fMinEventPlane = ep - TMath::Pi() / 4;
631 fMaxEventPlane = ep + TMath::Pi() / 4;
634 AliWarning("Could not set event plane limits because EMCal geometry was not loaded!");
638 //Load all requested track branches - each container knows name already
639 for (Int_t i =0; i<fParticleCollArray.GetEntriesFast(); i++) {
640 AliParticleContainer *cont = static_cast<AliParticleContainer*>(fParticleCollArray.At(i));
641 cont->SetArray(InputEvent());
644 if (fParticleCollArray.GetEntriesFast()>0) {
645 fTracks = GetParticleArray(0);
647 AliError(Form("%s: Could not retrieve first track branch!", GetName()));
652 //Load all requested cluster branches - each container knows name already
653 for (Int_t i =0; i<fClusterCollArray.GetEntriesFast(); i++) {
654 AliClusterContainer *cont = static_cast<AliClusterContainer*>(fClusterCollArray.At(i));
655 cont->SetArray(InputEvent());
658 if (fClusterCollArray.GetEntriesFast()>0) {
659 fCaloClusters = GetClusterArray(0);
660 if (!fCaloClusters) {
661 AliError(Form("%s: Could not retrieve first cluster branch!", GetName()));
666 if (!fCaloCellsName.IsNull() && !fCaloCells) {
667 fCaloCells = dynamic_cast<AliVCaloCells*>(InputEvent()->FindListObject(fCaloCellsName));
669 AliError(Form("%s: Could not retrieve cells %s!", GetName(), fCaloCellsName.Data()));
674 if (!fCaloTriggersName.IsNull() && !fCaloTriggers) {
675 fCaloTriggers = dynamic_cast<AliVCaloTrigger*>(InputEvent()->FindListObject(fCaloTriggersName));
676 if (!fCaloTriggers) {
677 AliError(Form("%s: Could not retrieve calo triggers %s!", GetName(), fCaloTriggersName.Data()));
682 if (!fCaloTriggerPatchInfoName.IsNull() && !fTriggerPatchInfo) {
683 fTriggerPatchInfo = GetArrayFromEvent(fCaloTriggerPatchInfoName.Data(),"AliEmcalTriggerPatchInfo");
684 if (!fTriggerPatchInfo) {
685 AliError(Form("%s: Could not retrieve calo trigger patch info %s!", GetName(), fCaloTriggerPatchInfoName.Data()));
691 fInitialized = kTRUE;
694 //_____________________________________________________
695 AliAnalysisTaskEmcal::BeamType AliAnalysisTaskEmcal::GetBeamType()
697 // Get beam type : pp-AA-pA
698 // ESDs have it directly, AODs get it from hardcoded run number ranges
700 if (fForceBeamType != kNA)
701 return fForceBeamType;
703 AliESDEvent *esd = dynamic_cast<AliESDEvent*>(InputEvent());
705 const AliESDRun *run = esd->GetESDRun();
706 TString beamType = run->GetBeamType();
707 if (beamType == "p-p")
709 else if (beamType == "A-A")
711 else if (beamType == "p-A")
716 Int_t runNumber = InputEvent()->GetRunNumber();
717 if ((runNumber >= 136851 && runNumber <= 139517) || // LHC10h
718 (runNumber >= 166529 && runNumber <= 170593)) { // LHC11h
720 } else if ((runNumber>=188365 && runNumber <= 188366) || // LHC12g
721 (runNumber >= 195344 && runNumber <= 196608)) { // LHC13b-f
729 //________________________________________________________________________
730 ULong_t AliAnalysisTaskEmcal::GetTriggerList()
732 if (!fTriggerPatchInfo)
735 //number of patches in event
736 Int_t nPatch = fTriggerPatchInfo->GetEntries();
738 //loop over patches to define trigger type of event
743 AliEmcalTriggerPatchInfo *patch;
744 for (Int_t iPatch = 0; iPatch < nPatch; iPatch++) {
745 patch = (AliEmcalTriggerPatchInfo*)fTriggerPatchInfo->At( iPatch );
746 if (patch->IsGammaHigh()) nG1++;
747 if (patch->IsGammaLow()) nG2++;
748 if (patch->IsJetHigh()) nJ1++;
749 if (patch->IsJetLow()) nJ2++;
752 AliDebug(2, "Patch summary: ");
753 AliDebug(2, Form("Number of patches: %d", nPatch));
754 AliDebug(2, Form("Jet: low[%d], high[%d]" ,nJ2, nJ1));
755 AliDebug(2, Form("Gamma: low[%d], high[%d]" ,nG2, nG1));
769 //________________________________________________________________________
770 Bool_t AliAnalysisTaskEmcal::HasTriggerType(Int_t trigger)
772 // Check if event has a given trigger type
774 return fTriggers == 0;
776 return (fTriggers & trigger);
779 //________________________________________________________________________
780 Bool_t AliAnalysisTaskEmcal::IsEventSelected()
782 // Check if event is selected
784 if (fOffTrigger != AliVEvent::kAny) {
786 const AliESDEvent *eev = dynamic_cast<const AliESDEvent*>(InputEvent());
788 res = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
790 const AliAODEvent *aev = dynamic_cast<const AliAODEvent*>(InputEvent());
792 res = ((AliVAODHeader*)aev->GetHeader())->GetOfflineTrigger();
795 if ((res & fOffTrigger) == 0) {
796 if (fGeneralHistograms) fHistEventRejection->Fill("PhysSel",1);
801 if (!fTrigClass.IsNull()) {
803 const AliESDEvent *eev = dynamic_cast<const AliESDEvent*>(InputEvent());
805 fired = eev->GetFiredTriggerClasses();
807 const AliAODEvent *aev = dynamic_cast<const AliAODEvent*>(InputEvent());
809 fired = aev->GetFiredTriggerClasses();
812 if (!fired.Contains("-B-")) {
813 if (fGeneralHistograms) fHistEventRejection->Fill("trigger",1);
816 TObjArray *arr = fTrigClass.Tokenize("|");
818 if (fGeneralHistograms) fHistEventRejection->Fill("trigger",1);
822 for (Int_t i=0;i<arr->GetEntriesFast();++i) {
823 TObject *obj = arr->At(i);
826 if (fired.Contains(obj->GetName())) {
833 if (fGeneralHistograms) fHistEventRejection->Fill("trigger",1);
838 if (!(fTriggerTypeSel & kND)) {
839 if (!HasTriggerType(fTriggerTypeSel)) {
840 if (fGeneralHistograms) fHistEventRejection->Fill("trigTypeSel",1);
845 if ((fMinCent != -999) && (fMaxCent != -999)) {
846 if (fCent<fMinCent || fCent>fMaxCent) {
847 if (fGeneralHistograms) fHistEventRejection->Fill("Cent",1);
852 if (fUseAliAnaUtils) {
853 if (!fAliAnalysisUtils)
854 fAliAnalysisUtils = new AliAnalysisUtils();
855 fAliAnalysisUtils->SetMinVtxContr(2);
856 fAliAnalysisUtils->SetMaxVtxZ(999);
857 if(fMinVz<-10.) fMinVz = -10.;
858 if(fMinVz>10.) fMaxVz = 10.;
860 if (!fAliAnalysisUtils->IsVertexSelected2013pA(InputEvent())) {
861 if (fGeneralHistograms) fHistEventRejection->Fill("VtxSel2013pA",1);
865 if (fRejectPileup && fAliAnalysisUtils->IsPileUpEvent(InputEvent())) {
866 if (fGeneralHistograms) fHistEventRejection->Fill("PileUp",1);
870 if(fTklVsClusSPDCut && fAliAnalysisUtils->IsSPDClusterVsTrackletBG(InputEvent())) {
871 if (fGeneralHistograms) fHistEventRejection->Fill("Bkg evt",1);
876 if ((fMinVz != -999) && (fMaxVz != -999)) {
877 if (fNVertCont == 0 ) {
878 if (fGeneralHistograms) fHistEventRejection->Fill("vertex contr.",1);
881 Double_t vz = fVertex[2];
882 if (vz<fMinVz || vz>fMaxVz) {
883 if (fGeneralHistograms) fHistEventRejection->Fill("Vz",1);
888 if (fMinPtTrackInEmcal > 0 && fGeom) {
889 Bool_t trackInEmcalOk = kFALSE;
890 Int_t ntracks = GetNParticles(0);
891 for (Int_t i = 0; i < ntracks; i++) {
892 AliVParticle *track = GetAcceptParticleFromArray(i,0);
896 Double_t phiMin = fGeom->GetArm1PhiMin() * TMath::DegToRad();
897 Double_t phiMax = fGeom->GetArm1PhiMax() * TMath::DegToRad();
898 Int_t runNumber = InputEvent()->GetRunNumber();
899 if (runNumber>=177295 && runNumber<=197470) { //small SM masked in 2012 and 2013
901 phiMax = TMath::Pi();
904 if (track->Eta() < fGeom->GetArm1EtaMin() || track->Eta() > fGeom->GetArm1EtaMax() || track->Phi() < phiMin || track->Phi() > phiMax)
906 if (track->Pt() > fMinPtTrackInEmcal) {
907 trackInEmcalOk = kTRUE;
911 if (!trackInEmcalOk) {
912 if (fGeneralHistograms) fHistEventRejection->Fill("trackInEmcal",1);
917 if (fMinNTrack > 0) {
918 Int_t nTracksAcc = 0;
919 Int_t ntracks = GetNParticles(0);
920 for (Int_t i = 0; i < ntracks; i++) {
921 AliVParticle *track = GetAcceptParticleFromArray(i,0);
924 if (track->Pt() > fTrackPtCut) {
926 if (nTracksAcc>=fMinNTrack)
930 if (nTracksAcc<fMinNTrack) {
931 if (fGeneralHistograms) fHistEventRejection->Fill("minNTrack",1);
936 if (!(fEPV0 > fMinEventPlane && fEPV0 <= fMaxEventPlane) &&
937 !(fEPV0 + TMath::Pi() > fMinEventPlane && fEPV0 + TMath::Pi() <= fMaxEventPlane) &&
938 !(fEPV0 - TMath::Pi() > fMinEventPlane && fEPV0 - TMath::Pi() <= fMaxEventPlane))
940 if (fGeneralHistograms) fHistEventRejection->Fill("EvtPlane",1);
944 if (fSelectPtHardBin != -999 && fSelectPtHardBin != fPtHardBin) {
945 if (fGeneralHistograms) fHistEventRejection->Fill("SelPtHardBin",1);
952 //________________________________________________________________________
953 TClonesArray *AliAnalysisTaskEmcal::GetArrayFromEvent(const char *name, const char *clname)
955 // Get array from event.
957 TClonesArray *arr = 0;
959 if (!sname.IsNull()) {
960 arr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(sname));
962 AliWarning(Form("%s: Could not retrieve array with name %s!", GetName(), name));
972 TString objname(arr->GetClass()->GetName());
974 if (!cls.InheritsFrom(clname)) {
975 AliWarning(Form("%s: Objects of type %s in %s are not inherited from %s!",
976 GetName(), cls.GetName(), name, clname));
982 //________________________________________________________________________
983 Bool_t AliAnalysisTaskEmcal::RetrieveEventObjects()
985 // Retrieve objects from event.
992 const AliVVertex *vert = InputEvent()->GetPrimaryVertex();
994 vert->GetXYZ(fVertex);
995 fNVertCont = vert->GetNContributors();
998 fBeamType = GetBeamType();
1000 if (fBeamType == kAA || fBeamType == kpA ) {
1001 AliCentrality *aliCent = InputEvent()->GetCentrality();
1003 fCent = aliCent->GetCentralityPercentile(fCentEst.Data());
1004 if (fNcentBins==4) {
1005 if (fCent >= 0 && fCent < 10) fCentBin = 0;
1006 else if (fCent >= 10 && fCent < 30) fCentBin = 1;
1007 else if (fCent >= 30 && fCent < 50) fCentBin = 2;
1008 else if (fCent >= 50 && fCent <= 100) fCentBin = 3;
1010 AliWarning(Form("%s: Negative centrality: %f. Assuming 99", GetName(), fCent));
1011 fCentBin = fNcentBins-1;
1014 Double_t centWidth = (fMaxCent-fMinCent)/(Double_t)fNcentBins;
1016 fCentBin = TMath::FloorNint(fCent/centWidth);
1019 if (fCentBin>=fNcentBins) {
1020 AliWarning(Form("%s: fCentBin too large: cent = %f fCentBin = %d. Assuming 99", GetName(),fCent,fCentBin));
1021 fCentBin = fNcentBins-1;
1025 AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
1026 fCentBin = fNcentBins-1;
1028 AliEventplane *aliEP = InputEvent()->GetEventplane();
1030 fEPV0 = aliEP->GetEventplane("V0" ,InputEvent());
1031 fEPV0A = aliEP->GetEventplane("V0A",InputEvent());
1032 fEPV0C = aliEP->GetEventplane("V0C",InputEvent());
1034 AliWarning(Form("%s: Could not retrieve event plane information!", GetName()));
1044 fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
1045 if (!fPythiaHeader) {
1047 AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
1050 for (UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++) {
1051 fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
1052 if (fPythiaHeader) break;
1058 if (fPythiaHeader) {
1059 fPtHard = fPythiaHeader->GetPtHard();
1061 const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
1062 const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
1063 for (fPtHardBin = 0; fPtHardBin < 11; fPtHardBin++) {
1064 if (fPtHard >= ptHardLo[fPtHardBin] && fPtHard < ptHardHi[fPtHardBin])
1068 fXsection = fPythiaHeader->GetXsection();
1069 fNTrials = fPythiaHeader->Trials();
1073 fTriggers = GetTriggerList();
1074 fMainTriggerPatchSet = kFALSE;
1075 memset(fMainTriggerPatch, 0, sizeof(AliEmcalTriggerPatchInfo *) * (kNType - 1));
1080 //________________________________________________________________________
1081 AliParticleContainer* AliAnalysisTaskEmcal::AddParticleContainer(const char *n)
1083 // Add particle container
1084 // will be called in AddTask macro
1086 TString tmp = TString(n);
1087 if (tmp.IsNull()) return 0;
1089 AliParticleContainer *cont = 0x0;
1090 cont = new AliParticleContainer();
1091 cont->SetArrayName(n);
1092 TString contName = cont->GetArrayName();
1094 fParticleCollArray.Add(cont);
1099 //________________________________________________________________________
1100 AliClusterContainer* AliAnalysisTaskEmcal::AddClusterContainer(const char *n)
1102 // Add cluster container
1103 // will be called in AddTask macro
1105 TString tmp = TString(n);
1106 if (tmp.IsNull()) return 0;
1108 AliClusterContainer *cont = 0x0;
1109 cont = new AliClusterContainer();
1110 cont->SetArrayName(n);
1112 fClusterCollArray.Add(cont);
1117 //________________________________________________________________________
1118 AliParticleContainer* AliAnalysisTaskEmcal::GetParticleContainer(Int_t i) const
1120 // Get i^th particle container
1122 if (i<0 || i>fParticleCollArray.GetEntriesFast()) return 0;
1123 AliParticleContainer *cont = static_cast<AliParticleContainer*>(fParticleCollArray.At(i));
1127 //________________________________________________________________________
1128 AliClusterContainer* AliAnalysisTaskEmcal::GetClusterContainer(Int_t i) const
1130 // Get i^th cluster container
1132 if (i<0 || i>fClusterCollArray.GetEntriesFast()) return 0;
1133 AliClusterContainer *cont = static_cast<AliClusterContainer*>(fClusterCollArray.At(i));
1137 //________________________________________________________________________
1138 AliParticleContainer* AliAnalysisTaskEmcal::GetParticleContainer(const char *name) const
1140 // Get particle container with name
1142 AliParticleContainer *cont = static_cast<AliParticleContainer*>(fParticleCollArray.FindObject(name));
1146 //________________________________________________________________________
1147 AliClusterContainer* AliAnalysisTaskEmcal::GetClusterContainer(const char *name) const
1149 // Get cluster container with name
1151 AliClusterContainer *cont = static_cast<AliClusterContainer*>(fClusterCollArray.FindObject(name));
1155 //________________________________________________________________________
1156 TClonesArray* AliAnalysisTaskEmcal::GetParticleArray(Int_t i) const
1158 // Get i^th TClonesArray with AliVParticle
1160 AliParticleContainer *cont = GetParticleContainer(i);
1162 AliError(Form("%s: Particle container %d not found",GetName(),i));
1165 TString contName = cont->GetArrayName();
1166 return cont->GetArray();
1169 //________________________________________________________________________
1170 TClonesArray* AliAnalysisTaskEmcal::GetClusterArray(Int_t i) const
1172 // Get i^th TClonesArray with AliVCluster
1174 AliClusterContainer *cont = GetClusterContainer(i);
1176 AliError(Form("%s:Cluster container %d not found",GetName(),i));
1179 return cont->GetArray();
1182 //________________________________________________________________________
1183 AliVParticle* AliAnalysisTaskEmcal::GetAcceptParticleFromArray(Int_t p, Int_t c) const
1185 // Get particle p if accepted from container c
1186 // If particle not accepted return 0
1188 AliParticleContainer *cont = GetParticleContainer(c);
1190 AliError(Form("%s: Particle container %d not found",GetName(),c));
1193 AliVParticle *vp = cont->GetAcceptParticle(p);
1198 //________________________________________________________________________
1199 AliVCluster* AliAnalysisTaskEmcal::GetAcceptClusterFromArray(Int_t cl, Int_t c) const
1201 // Get particle p if accepted from container c
1202 // If particle not accepted return 0
1204 AliClusterContainer *cont = GetClusterContainer(c);
1206 AliError(Form("%s: Cluster container %d not found",GetName(),c));
1209 AliVCluster *vc = cont->GetAcceptCluster(cl);
1214 //________________________________________________________________________
1215 Int_t AliAnalysisTaskEmcal::GetNParticles(Int_t i) const
1217 // Get number of entries in particle array i
1219 AliParticleContainer *cont = GetParticleContainer(i);
1221 AliError(Form("%s: Particle container %d not found",GetName(),i));
1224 return cont->GetNEntries();
1227 //________________________________________________________________________
1228 Int_t AliAnalysisTaskEmcal::GetNClusters(Int_t i) const
1230 // Get number of entries in cluster array i
1232 AliClusterContainer *cont = GetClusterContainer(i);
1234 AliError(Form("%s: Cluster container %d not found",GetName(),i));
1237 return cont->GetNEntries();
1240 //________________________________________________________________________
1241 AliEmcalTriggerPatchInfo* AliAnalysisTaskEmcal::GetMainTriggerPatch(Int_t types)
1243 //get main trigger match; if not known yet, look for it and cache
1245 // In order to be able to select more than one class at the same time (i.e. low and
1246 // high threshold trigger), the argument must be a bit string, however using the enum
1247 // of trigger classes in order to define the selection
1249 // Get from cache (if already set)
1250 if(fMainTriggerPatchSet){
1251 return ApplyMainTriggerSelection(types);
1254 if (!fTriggerPatchInfo) {
1255 AliError(Form("%s: fTriggerPatchInfo not available",GetName()));
1259 //number of patches in event
1260 Int_t nPatch = fTriggerPatchInfo->GetEntries();
1262 //extract main trigger patch(es)
1263 AliEmcalTriggerPatchInfo *patch;
1264 for (Int_t iPatch = 0; iPatch < nPatch; iPatch++) {
1266 patch = (AliEmcalTriggerPatchInfo*)fTriggerPatchInfo->At( iPatch );
1267 if (patch->IsMainTrigger()) {
1268 if(patch->IsJetHigh()) fMainTriggerPatch[kPosJ1] = patch;
1269 if(patch->IsJetLow()) fMainTriggerPatch[kPosJ2] = patch;
1270 if(patch->IsGammaHigh()) fMainTriggerPatch[kPosG1] = patch;
1271 if(patch->IsGammaLow()) fMainTriggerPatch[kPosG2] = patch;
1272 if(patch->IsLevel0()) fMainTriggerPatch[kPosL0] = patch;
1275 fMainTriggerPatchSet = kTRUE;
1277 return ApplyMainTriggerSelection(types);
1280 AliEmcalTriggerPatchInfo *AliAnalysisTaskEmcal::ApplyMainTriggerSelection(Int_t selectionBitmap) const {
1282 // Select main trigger patch for a given combination of trigger classes
1283 // If the bitmap contains only one class, the selection is trivial
1284 // If the bitmap contains many classes, the selection takes the highest of the patchess according to the
1287 AliEmcalTriggerPatchInfo *highest(NULL);
1288 if(selectionBitmap & kJ1){
1289 highest = fMainTriggerPatch[kPosJ1];
1291 if(selectionBitmap & kJ2){
1292 if(!highest) highest = fMainTriggerPatch[kPosJ2];
1294 if(fMainTriggerPatch[kPosJ2] && fMainTriggerPatch[kPosJ2]->GetADCAmp() > highest->GetADCAmp())
1295 highest = fMainTriggerPatch[kPosJ2];
1299 if(selectionBitmap & kG1){
1300 if(!highest) highest = fMainTriggerPatch[kPosG1];
1302 if(fMainTriggerPatch[kPosG1] && fMainTriggerPatch[kPosG1]->GetADCAmp() > highest->GetADCAmp())
1303 highest = fMainTriggerPatch[kPosG1];
1307 if(selectionBitmap & kG2){
1308 if(!highest) highest = fMainTriggerPatch[kPosG2];
1310 if(fMainTriggerPatch[kPosG2] && fMainTriggerPatch[kPosG2]->GetADCAmp() > highest->GetADCAmp())
1311 highest = fMainTriggerPatch[kPosG2];
1314 if(selectionBitmap & kL0){
1315 if(!highest) highest = fMainTriggerPatch[kPosL0];
1317 if(fMainTriggerPatch[kPosL0] && fMainTriggerPatch[kPosL0]->GetADCAmp() > highest->GetADCAmp())
1318 highest = fMainTriggerPatch[kPosL0];
1325 //________________________________________________________________________
1326 void AliAnalysisTaskEmcal::AddObjectToEvent(TObject *obj)
1328 // Add object to event
1330 if (!(InputEvent()->FindListObject(obj->GetName()))) {
1331 InputEvent()->AddObject(obj);
1333 AliFatal(Form("%s: Container with name %s already present. Aborting", GetName(), obj->GetName()));
1337 //________________________________________________________________________
1338 void AliAnalysisTaskEmcal::GenerateFixedBinArray(Int_t n, Double_t min, Double_t max, Double_t* array) const
1340 Double_t binWidth = (max-min)/n;
1342 for (Int_t i = 1; i <= n; i++) {
1343 array[i] = array[i-1]+binWidth;
1347 //________________________________________________________________________
1348 Double_t* AliAnalysisTaskEmcal::GenerateFixedBinArray(Int_t n, Double_t min, Double_t max) const
1350 Double_t *array = new Double_t[n+1];
1352 GenerateFixedBinArray(n, min, max, array);
1357 //________________________________________________________________________
1358 void AliAnalysisTaskEmcal::SetRejectionReasonLabels(TAxis* axis)
1360 axis->SetBinLabel(1, "NullObject");
1361 axis->SetBinLabel(2, "Pt");
1362 axis->SetBinLabel(3, "Acceptance");
1363 axis->SetBinLabel(4, "BitMap");
1364 axis->SetBinLabel(5, "Bit4");
1365 axis->SetBinLabel(6, "Bit5");
1366 axis->SetBinLabel(7, "Bit6");
1367 axis->SetBinLabel(8, "Bit7");
1368 axis->SetBinLabel(9, "MCFlag");
1369 axis->SetBinLabel(10, "MCGenerator");
1370 axis->SetBinLabel(11, "ChargeCut");
1371 axis->SetBinLabel(12, "Bit11");
1372 axis->SetBinLabel(13, "Bit12");
1373 axis->SetBinLabel(14, "IsEMCal");
1374 axis->SetBinLabel(15, "Time");
1375 axis->SetBinLabel(16, "Energy");
1376 axis->SetBinLabel(17, "Bit16");
1377 axis->SetBinLabel(18, "Bit17");
1378 axis->SetBinLabel(19, "Area");
1379 axis->SetBinLabel(20, "AreaEmc");
1380 axis->SetBinLabel(21, "ZLeadingCh");
1381 axis->SetBinLabel(22, "ZLeadingEmc");
1382 axis->SetBinLabel(23, "NEF");
1383 axis->SetBinLabel(24, "MinLeadPt");
1384 axis->SetBinLabel(25, "MaxTrackPt");
1385 axis->SetBinLabel(26, "MaxClusterPt");
1386 axis->SetBinLabel(27, "Flavour");
1387 axis->SetBinLabel(28, "TagStatus");
1388 axis->SetBinLabel(29, "Bit28");
1389 axis->SetBinLabel(30, "Bit29");
1390 axis->SetBinLabel(31, "Bit30");
1391 axis->SetBinLabel(32, "Bit31");