]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG/EMCAL/AliAnalysisTaskEmcal.cxx
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWG / EMCAL / AliAnalysisTaskEmcal.cxx
index bcf88b4d161697dd22ad78db0084ed93f01dc784..aeb6b7817165ba3be33ba42303a6c333677caa45 100644 (file)
@@ -1,4 +1,3 @@
-// $Id$
 //
 // Emcal base analysis task.
 //
@@ -58,6 +57,8 @@ AliAnalysisTaskEmcal::AliAnalysisTaskEmcal() :
   fTrackPtCut(0),
   fMinNTrack(0),
   fUseAliAnaUtils(kFALSE),
+  fRejectPileup(kFALSE),
+  fTklVsClusSPDCut(kFALSE),
   fAliAnalysisUtils(0x0),
   fOffTrigger(AliVEvent::kAny),
   fTrigClass(),
@@ -76,6 +77,8 @@ AliAnalysisTaskEmcal::AliAnalysisTaskEmcal() :
   fMinMCLabel(0),
   fMCLabelShift(0),
   fNcentBins(4),
+  fNeedEmcalGeom(kTRUE),
+  fIsEsd(kFALSE),
   fGeom(0),
   fTracks(0),
   fCaloClusters(0),
@@ -93,16 +96,18 @@ AliAnalysisTaskEmcal::AliAnalysisTaskEmcal() :
   fPtHard(0),
   fPtHardBin(0),
   fNTrials(0),
+  fXsection(0),
   fParticleCollArray(),
   fClusterCollArray(),
-  fMainTriggerPatch(0x0),
-  fTriggerType(kND),
+  fTriggers(0),
   fOutput(0),
+  fHistEventCount(0),
   fHistTrialsAfterSel(0),
   fHistEventsAfterSel(0),
+  fHistXsectionAfterSel(0),
   fHistTrials(0),
-  fHistXsection(0),
   fHistEvents(0),
+  fHistXsection(0),
   fHistPtHard(0),
   fHistCentrality(0),
   fHistZVertex(0),
@@ -136,6 +141,8 @@ AliAnalysisTaskEmcal::AliAnalysisTaskEmcal(const char *name, Bool_t histo) :
   fTrackPtCut(0),
   fMinNTrack(0),
   fUseAliAnaUtils(kFALSE),
+  fRejectPileup(kFALSE),
+  fTklVsClusSPDCut(kFALSE),
   fAliAnalysisUtils(0x0),
   fOffTrigger(AliVEvent::kAny),
   fTrigClass(),
@@ -154,6 +161,8 @@ AliAnalysisTaskEmcal::AliAnalysisTaskEmcal(const char *name, Bool_t histo) :
   fMinMCLabel(0),
   fMCLabelShift(0),
   fNcentBins(4),
+  fNeedEmcalGeom(kTRUE),
+  fIsEsd(kFALSE),
   fGeom(0),
   fTracks(0),
   fCaloClusters(0),
@@ -171,16 +180,18 @@ AliAnalysisTaskEmcal::AliAnalysisTaskEmcal(const char *name, Bool_t histo) :
   fPtHard(0),
   fPtHardBin(0),
   fNTrials(0),
+  fXsection(0),
   fParticleCollArray(),
   fClusterCollArray(),
-  fMainTriggerPatch(0x0),
-  fTriggerType(kND),
+  fTriggers(0),
   fOutput(0),
+  fHistEventCount(0),
   fHistTrialsAfterSel(0),
   fHistEventsAfterSel(0),
+  fHistXsectionAfterSel(0),
   fHistTrials(0),
-  fHistXsection(0),
   fHistEvents(0),
+  fHistXsection(0),
   fHistPtHard(0),
   fHistCentrality(0),
   fHistZVertex(0),
@@ -253,6 +264,26 @@ void AliAnalysisTaskEmcal::SetTrackPhiLimits(Double_t min, Double_t max, Int_t c
 void AliAnalysisTaskEmcal::UserCreateOutputObjects()
 {
   // Create user output.
+
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (mgr) {
+    AliVEventHandler *evhand = mgr->GetInputEventHandler();
+    if (evhand) {
+      if (evhand->InheritsFrom("AliESDInputHandler")) {
+        fIsEsd = kTRUE;
+      }
+      else {
+        fIsEsd = kFALSE;        
+      }
+    }
+    else {
+      AliError("Event handler not found!");
+    }
+  }
+  else {
+    AliError("Analysis manager not found!");
+  }  
+
   if (!fCreateHisto)
     return;
 
@@ -276,22 +307,27 @@ void AliAnalysisTaskEmcal::UserCreateOutputObjects()
     fHistEventsAfterSel->GetXaxis()->SetTitle("p_{T} hard bin");
     fHistEventsAfterSel->GetYaxis()->SetTitle("total events");
     fOutput->Add(fHistEventsAfterSel);
+
+    fHistXsectionAfterSel = new TProfile("fHistXsectionAfterSel", "fHistXsectionAfterSel", 11, 0, 11);
+    fHistXsectionAfterSel->GetXaxis()->SetTitle("p_{T} hard bin");
+    fHistXsectionAfterSel->GetYaxis()->SetTitle("xsection");
+    fOutput->Add(fHistXsectionAfterSel);
     
     fHistTrials = new TH1F("fHistTrials", "fHistTrials", 11, 0, 11);
     fHistTrials->GetXaxis()->SetTitle("p_{T} hard bin");
     fHistTrials->GetYaxis()->SetTitle("trials");
     fOutput->Add(fHistTrials);
 
-    fHistXsection = new TProfile("fHistXsection", "fHistXsection", 11, 0, 11);
-    fHistXsection->GetXaxis()->SetTitle("p_{T} hard bin");
-    fHistXsection->GetYaxis()->SetTitle("xsection");
-    fOutput->Add(fHistXsection);
-
     fHistEvents = new TH1F("fHistEvents", "fHistEvents", 11, 0, 11);
     fHistEvents->GetXaxis()->SetTitle("p_{T} hard bin");
     fHistEvents->GetYaxis()->SetTitle("total events");
     fOutput->Add(fHistEvents);
 
+    fHistXsection = new TProfile("fHistXsection", "fHistXsection", 11, 0, 11);
+    fHistXsection->GetXaxis()->SetTitle("p_{T} hard bin");
+    fHistXsection->GetYaxis()->SetTitle("xsection");
+    fOutput->Add(fHistXsection);
+
     const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
     const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
     
@@ -325,22 +361,29 @@ void AliAnalysisTaskEmcal::UserCreateOutputObjects()
   fHistEventPlane->GetYaxis()->SetTitle("counts");
   fOutput->Add(fHistEventPlane);
 
-  fHistEventRejection = new TH1I("fHistEventRejection","Reasons to reject event",20,0,20);
-  fHistEventRejection->Fill("PhysSel",0);
-  fHistEventRejection->Fill("trigger",0);
-  fHistEventRejection->Fill("trigTypeSel",0);
-  fHistEventRejection->Fill("Cent",0);
-  fHistEventRejection->Fill("vertex contr.",0);
-  fHistEventRejection->Fill("Vz",0);
-  fHistEventRejection->Fill("trackInEmcal",0);
-  fHistEventRejection->Fill("minNTrack",0);
-  fHistEventRejection->Fill("VtxSel2013pA",0);
-  fHistEventRejection->Fill("PileUp",0);
-  fHistEventRejection->Fill("EvtPlane",0);
-  fHistEventRejection->Fill("SelPtHardBin",0);
+  fHistEventRejection = new TH1F("fHistEventRejection","Reasons to reject event",20,0,20);
+  fHistEventRejection->GetXaxis()->SetBinLabel(1,"PhysSel");
+  fHistEventRejection->GetXaxis()->SetBinLabel(2,"trigger");
+  fHistEventRejection->GetXaxis()->SetBinLabel(3,"trigTypeSel");
+  fHistEventRejection->GetXaxis()->SetBinLabel(4,"Cent");
+  fHistEventRejection->GetXaxis()->SetBinLabel(5,"vertex contr.");
+  fHistEventRejection->GetXaxis()->SetBinLabel(6,"Vz");
+  fHistEventRejection->GetXaxis()->SetBinLabel(7,"trackInEmcal");
+  fHistEventRejection->GetXaxis()->SetBinLabel(8,"minNTrack");
+  fHistEventRejection->GetXaxis()->SetBinLabel(9,"VtxSel2013pA");
+  fHistEventRejection->GetXaxis()->SetBinLabel(10,"PileUp");
+  fHistEventRejection->GetXaxis()->SetBinLabel(11,"EvtPlane");
+  fHistEventRejection->GetXaxis()->SetBinLabel(12,"SelPtHardBin");
+  fHistEventRejection->GetXaxis()->SetBinLabel(13,"Bkg evt");
   fHistEventRejection->GetYaxis()->SetTitle("counts");
   fOutput->Add(fHistEventRejection);
 
+  fHistEventCount = new TH1F("fHistEventCount","fHistEventCount",2,0,2);
+  fHistEventCount->GetXaxis()->SetBinLabel(1,"Accepted");
+  fHistEventCount->GetXaxis()->SetBinLabel(2,"Rejected");
+  fHistEventCount->GetYaxis()->SetTitle("counts");
+  fOutput->Add(fHistEventCount);
+
   PostData(1, fOutput);
 }
 
@@ -348,8 +391,9 @@ void AliAnalysisTaskEmcal::UserCreateOutputObjects()
 Bool_t AliAnalysisTaskEmcal::FillGeneralHistograms()
 {
   if (fIsPythia) {
-    fHistEventsAfterSel->SetBinContent(fPtHardBin + 1, fHistEventsAfterSel->GetBinContent(fPtHardBin + 1) + 1);
-    fHistTrialsAfterSel->SetBinContent(fPtHardBin + 1, fHistTrialsAfterSel->GetBinContent(fPtHardBin + 1) + fNTrials);
+    fHistEventsAfterSel->Fill(fPtHardBin, 1);
+    fHistTrialsAfterSel->Fill(fPtHardBin, fNTrials);
+    fHistXsectionAfterSel->Fill(fPtHardBin, fXsection);
     fHistPtHard->Fill(fPtHard);
   }
 
@@ -365,8 +409,6 @@ void AliAnalysisTaskEmcal::UserExec(Option_t *)
 {
   // Main loop, called for each event.
 
-  fMainTriggerPatch = NULL;
-
   if (!fInitialized)
     ExecOnce();
 
@@ -376,8 +418,13 @@ void AliAnalysisTaskEmcal::UserExec(Option_t *)
   if (!RetrieveEventObjects())
     return;
 
-  if (!IsEventSelected()) 
+  if (IsEventSelected()) {
+    if (fGeneralHistograms) fHistEventCount->Fill("Accepted",1);
+  }
+  else {
+    if (fGeneralHistograms) fHistEventCount->Fill("Rejected",1);
     return;
+  }
 
   if (fGeneralHistograms && fCreateHisto) {
     if (!FillGeneralHistograms())
@@ -413,7 +460,6 @@ Bool_t AliAnalysisTaskEmcal::AcceptCluster(AliVCluster *clus, Int_t c) const
   }
 
   return cont->AcceptCluster(clus);
-
 }
 
 //________________________________________________________________________
@@ -451,8 +497,7 @@ Bool_t AliAnalysisTaskEmcal::PythiaInfoFromFile(const char* currFile, Float_t &f
     Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
     Ssiz_t pos2 = file.Index(".root",5,TString::kExact);
     file.Replace(pos+1,pos2-pos1,"");
-  }
-  else {
+  } else {
     // not an archive take the basename....
     file.ReplaceAll(gSystem->BaseName(file.Data()),"");
   }
@@ -479,8 +524,7 @@ Bool_t AliAnalysisTaskEmcal::PythiaInfoFromFile(const char* currFile, Float_t &f
     if (!fxsec) {
        // not a severe condition but inciate that we have no information
       return kFALSE;
-    }
-    else{
+    } else {
       // find the tlist we want to be independtent of the name so use the Tkey
       TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0); 
       if (!key) {
@@ -496,8 +540,7 @@ Bool_t AliAnalysisTaskEmcal::PythiaInfoFromFile(const char* currFile, Float_t &f
       fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
       fxsec->Close();
     }
-  } // no tree pyxsec.root
-  else {
+  } else { // no tree pyxsec.root
     TTree *xtree = (TTree*)fxsec->Get("Xsection");
     if (!xtree) {
       fxsec->Close();
@@ -518,6 +561,8 @@ Bool_t AliAnalysisTaskEmcal::PythiaInfoFromFile(const char* currFile, Float_t &f
 //________________________________________________________________________
 Bool_t AliAnalysisTaskEmcal::UserNotify()
 {
+  // Called when file changes.
+
   if (!fIsPythia || !fGeneralHistograms || !fCreateHisto)
     return kTRUE;
 
@@ -527,9 +572,9 @@ Bool_t AliAnalysisTaskEmcal::UserNotify()
     return kFALSE;
   }
 
-  Float_t xsection = 0;
-  Float_t trials   = 0;
-  Int_t   pthard   = 0;
+  Float_t xsection    = 0;
+  Float_t trials      = 0;
+  Int_t   pthardbin   = 0;
 
   TFile *curfile = tree->GetCurrentFile();
   if (!curfile) {
@@ -538,19 +583,18 @@ Bool_t AliAnalysisTaskEmcal::UserNotify()
   }
 
   TChain *chain = dynamic_cast<TChain*>(tree);
-  if (chain)
-    tree = chain->GetTree();
+  if (chain) tree = chain->GetTree();
 
   Int_t nevents = tree->GetEntriesFast();
 
-  PythiaInfoFromFile(curfile->GetName(), xsection, trials, pthard);
+  PythiaInfoFromFile(curfile->GetName(), xsection, trials, pthardbin);
 
   // TODO: Workaround
-  if ((pthard < 0) || (pthard > 10))
-    pthard = 0;
-  fHistTrials->Fill(pthard, trials);
-  fHistXsection->Fill(pthard, xsection);
-  fHistEvents->Fill(pthard, nevents);
+  if ((pthardbin < 0) || (pthardbin > 10)) pthardbin = 0;
+
+  fHistTrials->Fill(pthardbin, trials);
+  fHistXsection->Fill(pthardbin, xsection);
+  fHistEvents->Fill(pthardbin, nevents);
 
   return kTRUE;
 }
@@ -565,16 +609,24 @@ void AliAnalysisTaskEmcal::ExecOnce()
     return;
   }
 
-  fGeom = AliEMCALGeometry::GetInstance();
-  if (!fGeom) {
-    AliError(Form("%s: Can not create geometry", GetName()));
-    return;
+  if (fNeedEmcalGeom) {
+    fGeom = AliEMCALGeometry::GetInstance();
+    if (!fGeom) {
+      AliError(Form("%s: Can not create geometry", GetName()));
+      return;
+    }
   }
 
+  
   if (fEventPlaneVsEmcal >= 0) {
-    Double_t ep = (fGeom->GetArm1PhiMax() + fGeom->GetArm1PhiMin()) / 2 * TMath::DegToRad() + fEventPlaneVsEmcal - TMath::Pi();
-    fMinEventPlane = ep - TMath::Pi() / 4;
-    fMaxEventPlane = ep + TMath::Pi() / 4;
+    if (fGeom) {
+      Double_t ep = (fGeom->GetArm1PhiMax() + fGeom->GetArm1PhiMin()) / 2 * TMath::DegToRad() + fEventPlaneVsEmcal - TMath::Pi();
+      fMinEventPlane = ep - TMath::Pi() / 4;
+      fMaxEventPlane = ep + TMath::Pi() / 4;
+    }
+    else {
+      AliWarning("Could not set event plane limits because EMCal geometry was not loaded!");
+    }
   }
 
   //Load all requested track branches - each container knows name already
@@ -657,13 +709,10 @@ AliAnalysisTaskEmcal::BeamType AliAnalysisTaskEmcal::GetBeamType()
   } else {
     Int_t runNumber = InputEvent()->GetRunNumber();
     if ((runNumber >= 136851 && runNumber <= 139517) ||  // LHC10h
-       (runNumber >= 166529 && runNumber <= 170593))    // LHC11h
-    { 
+       (runNumber >= 166529 && runNumber <= 170593)) {  // LHC11h
       return kAA;
-    } 
-    else if ((runNumber>=188365 && runNumber <= 188366) || // LHC12g
-            (runNumber >= 195344 && runNumber <= 196608))  // LHC13b-f
-    {
+    } else if ((runNumber>=188365 && runNumber <= 188366) ||   // LHC12g
+              (runNumber >= 195344 && runNumber <= 196608)) { // LHC13b-f
       return kpA;
     } else {
       return kpp;
@@ -671,34 +720,53 @@ AliAnalysisTaskEmcal::BeamType AliAnalysisTaskEmcal::GetBeamType()
   }  
 }
 
-//_____________________________________________________
-AliAnalysisTaskEmcal::TriggerType AliAnalysisTaskEmcal::GetTriggerType()
+//________________________________________________________________________
+ULong_t AliAnalysisTaskEmcal::GetTriggerList()
 {
-  // Get trigger type: kND, kJ1, kJ2
   if (!fTriggerPatchInfo)
-    return kND;
-  
+    return 0;
+
   //number of patches in event
   Int_t nPatch = fTriggerPatchInfo->GetEntries();
 
   //loop over patches to define trigger type of event
+  Int_t nG1 = 0;
+  Int_t nG2 = 0;
   Int_t nJ1 = 0;
   Int_t nJ2 = 0;
+  Int_t nL0 = 0;
   AliEmcalTriggerPatchInfo *patch;
   for (Int_t iPatch = 0; iPatch < nPatch; iPatch++) {
     patch = (AliEmcalTriggerPatchInfo*)fTriggerPatchInfo->At( iPatch );
+    if (patch->IsGammaHigh()) nG1++;
+    if (patch->IsGammaLow())  nG2++;
     if (patch->IsJetHigh()) nJ1++;
     if (patch->IsJetLow())  nJ2++;
+    if (patch->IsLevel0())  nL0++;
   }
 
-  if (nJ1>0) 
-    return kJ1;
-  else if (nJ2>0)
-    return kJ2;
-  else
-    return kND;
+  AliDebug(2, "Patch summary: ");
+  AliDebug(2, Form("Number of patches: %d", nPatch));
+  AliDebug(2, Form("Jet:   low[%d], high[%d]" ,nJ2, nJ1));
+  AliDebug(2, Form("Gamma: low[%d], high[%d]" ,nG2, nG1));
+
+  ULong_t triggers(0);
+  if (nL0>0) SETBIT(triggers, kL0);
+  if (nG1>0) SETBIT(triggers, kG1);
+  if (nG2>0) SETBIT(triggers, kG2);
+  if (nJ1>0) SETBIT(triggers, kJ1);
+  if (nJ2>0) SETBIT(triggers, kJ2);
+  return triggers;
+}
+
+//________________________________________________________________________
+Bool_t AliAnalysisTaskEmcal::HasTriggerType(TriggerType trigger)
+{
+  // Check if event has a given trigger type
+  if(trigger & kND){
+    return fTriggers == 0;
+  }
+  return TESTBIT(fTriggers, trigger);
 }
 
 //________________________________________________________________________
@@ -714,12 +782,11 @@ Bool_t AliAnalysisTaskEmcal::IsEventSelected()
     } else {
       const AliAODEvent *aev = dynamic_cast<const AliAODEvent*>(InputEvent());
       if (aev) {
-        res = aev->GetHeader()->GetOfflineTrigger();
+        res = ((AliVAODHeader*)aev->GetHeader())->GetOfflineTrigger();
       }
     }
     if ((res & fOffTrigger) == 0) {
-      if (fGeneralHistograms) 
-       fHistEventRejection->Fill("PhysSel",1);
+      if (fGeneralHistograms) fHistEventRejection->Fill("PhysSel",1);
       return kFALSE;
     }
   }
@@ -736,14 +803,12 @@ Bool_t AliAnalysisTaskEmcal::IsEventSelected()
       }
     }
     if (!fired.Contains("-B-")) {
-      if (fGeneralHistograms) 
-       fHistEventRejection->Fill("trigger",1);
+      if (fGeneralHistograms) fHistEventRejection->Fill("trigger",1);
       return kFALSE;
     }
     TObjArray *arr = fTrigClass.Tokenize("|");
     if (!arr) {
-      if (fGeneralHistograms) 
-       fHistEventRejection->Fill("trigger",1);
+      if (fGeneralHistograms) fHistEventRejection->Fill("trigger",1);
       return kFALSE;
     }
     Bool_t match = 0;
@@ -758,39 +823,57 @@ Bool_t AliAnalysisTaskEmcal::IsEventSelected()
     }
     delete arr;
     if (!match) {
-      if (fGeneralHistograms) 
-       fHistEventRejection->Fill("trigger",1);
+      if (fGeneralHistograms) fHistEventRejection->Fill("trigger",1);
       return kFALSE;
     }
   }
 
   if (fTriggerTypeSel != kND) {
-    if (fTriggerType != fTriggerTypeSel) {
-      if (fGeneralHistograms) 
-       fHistEventRejection->Fill("trigTypeSel",1);
+    if (!HasTriggerType(fTriggerTypeSel)) {
+      if (fGeneralHistograms) fHistEventRejection->Fill("trigTypeSel",1);
       return kFALSE;
     }
   }
-  
 
   if ((fMinCent != -999) && (fMaxCent != -999)) {
     if (fCent<fMinCent || fCent>fMaxCent) {
-      if (fGeneralHistograms) 
-       fHistEventRejection->Fill("Cent",1);
+      if (fGeneralHistograms) fHistEventRejection->Fill("Cent",1);
+      return kFALSE;
+    }
+  }
+
+  if (fUseAliAnaUtils) {
+    if (!fAliAnalysisUtils)
+      fAliAnalysisUtils = new AliAnalysisUtils();
+    fAliAnalysisUtils->SetMinVtxContr(2);
+    fAliAnalysisUtils->SetMaxVtxZ(999);
+    if(fMinVz<-10.) fMinVz = -10.; 
+    if(fMinVz>10.)  fMaxVz = 10.;
+
+    if (!fAliAnalysisUtils->IsVertexSelected2013pA(InputEvent())) {
+      if (fGeneralHistograms) fHistEventRejection->Fill("VtxSel2013pA",1);
+      return kFALSE;
+    }
+
+    if (fRejectPileup && fAliAnalysisUtils->IsPileUpEvent(InputEvent())) {
+      if (fGeneralHistograms) fHistEventRejection->Fill("PileUp",1);
+      return kFALSE;
+    }
+
+    if(fTklVsClusSPDCut && fAliAnalysisUtils->IsSPDClusterVsTrackletBG(InputEvent())) {
+      if (fGeneralHistograms) fHistEventRejection->Fill("Bkg evt",1);
       return kFALSE;
     }
   }
 
   if ((fMinVz != -999) && (fMaxVz != -999)) {
     if (fNVertCont == 0 ) {
-      if (fGeneralHistograms) 
-       fHistEventRejection->Fill("vertex contr.",1);
+      if (fGeneralHistograms) fHistEventRejection->Fill("vertex contr.",1);
       return kFALSE;
     }
     Double_t vz = fVertex[2];
     if (vz<fMinVz || vz>fMaxVz) {
-      if (fGeneralHistograms) 
-       fHistEventRejection->Fill("Vz",1);
+      if (fGeneralHistograms) fHistEventRejection->Fill("Vz",1);
       return kFALSE;
     }
   }
@@ -819,8 +902,7 @@ Bool_t AliAnalysisTaskEmcal::IsEventSelected()
       }
     }
     if (!trackInEmcalOk) {
-      if (fGeneralHistograms) 
-       fHistEventRejection->Fill("trackInEmcal",1);
+      if (fGeneralHistograms) fHistEventRejection->Fill("trackInEmcal",1);
       return kFALSE;
     }
   }
@@ -839,26 +921,7 @@ Bool_t AliAnalysisTaskEmcal::IsEventSelected()
       }
     }
     if (nTracksAcc<fMinNTrack) {
-      if (fGeneralHistograms) 
-       fHistEventRejection->Fill("minNTrack",1);
-      return kFALSE;
-    }
-  }
-
-  if (fUseAliAnaUtils) {
-    if (!fAliAnalysisUtils)
-      fAliAnalysisUtils = new AliAnalysisUtils();
-    fAliAnalysisUtils->SetMinVtxContr(2);
-    fAliAnalysisUtils->SetMaxVtxZ(10.);
-
-    if (!fAliAnalysisUtils->IsVertexSelected2013pA(InputEvent())) {
-      if (fGeneralHistograms) 
-       fHistEventRejection->Fill("VtxSel2013pA",1);
-      return kFALSE;
-    }
-
-    if (fAliAnalysisUtils->IsPileUpEvent(InputEvent())) {
-      fHistEventRejection->Fill("PileUp",1);
+      if (fGeneralHistograms) fHistEventRejection->Fill("minNTrack",1);
       return kFALSE;
     }
   }
@@ -867,17 +930,15 @@ Bool_t AliAnalysisTaskEmcal::IsEventSelected()
       !(fEPV0 + TMath::Pi() > fMinEventPlane && fEPV0 + TMath::Pi() <= fMaxEventPlane) &&
       !(fEPV0 - TMath::Pi() > fMinEventPlane && fEPV0 - TMath::Pi() <= fMaxEventPlane)) 
     {
-      if (fGeneralHistograms) 
-       fHistEventRejection->Fill("EvtPlane",1);
+      if (fGeneralHistograms) fHistEventRejection->Fill("EvtPlane",1);
       return kFALSE;
     }
 
   if (fSelectPtHardBin != -999 && fSelectPtHardBin != fPtHardBin)  {
-      if (fGeneralHistograms) 
-       fHistEventRejection->Fill("SelPtHardBin",1);
+      if (fGeneralHistograms) fHistEventRejection->Fill("SelPtHardBin",1);
       return kFALSE;
     }
-
+  
   return kTRUE;
 }
 
@@ -942,10 +1003,12 @@ Bool_t AliAnalysisTaskEmcal::RetrieveEventObjects()
          AliWarning(Form("%s: Negative centrality: %f. Assuming 99", GetName(), fCent));
          fCentBin = fNcentBins-1;
        }
-      }
-      else {
+      } else {
        Double_t centWidth = (fMaxCent-fMinCent)/(Double_t)fNcentBins;
-       fCentBin = TMath::FloorNint(fCent/centWidth);
+       if(centWidth>0.)
+         fCentBin = TMath::FloorNint(fCent/centWidth);
+       else 
+         fCentBin = 0;
        if (fCentBin>=fNcentBins) {
          AliWarning(Form("%s: fCentBin too large: cent = %f fCentBin = %d. Assuming 99", GetName(),fCent,fCentBin));
          fCentBin = fNcentBins-1;
@@ -953,7 +1016,7 @@ Bool_t AliAnalysisTaskEmcal::RetrieveEventObjects()
       }
     } else {
       AliWarning(Form("%s: Could not retrieve centrality information! Assuming 99", GetName()));
-      fCentBin = 3;
+      fCentBin = fNcentBins-1;
     }
     AliEventplane *aliEP = InputEvent()->GetEventplane();
     if (aliEP) {
@@ -987,26 +1050,27 @@ Bool_t AliAnalysisTaskEmcal::RetrieveEventObjects()
 
     if (fPythiaHeader) {
       fPtHard = fPythiaHeader->GetPtHard();
-    
+
       const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
       const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
       for (fPtHardBin = 0; fPtHardBin < 11; fPtHardBin++) {
        if (fPtHard >= ptHardLo[fPtHardBin] && fPtHard < ptHardHi[fPtHardBin])
          break;
       }
-    
+
+      fXsection = fPythiaHeader->GetXsection();
       fNTrials = fPythiaHeader->Trials();
     }
   }
 
-  fTriggerType = GetTriggerType();
+  fTriggers = GetTriggerList();
 
   return kTRUE;
 }
 
 //________________________________________________________________________
-AliParticleContainer* AliAnalysisTaskEmcal::AddParticleContainer(const char *n) {
-
+AliParticleContainer* AliAnalysisTaskEmcal::AddParticleContainer(const char *n) 
+{
   // Add particle container
   // will be called in AddTask macro
 
@@ -1024,8 +1088,8 @@ AliParticleContainer* AliAnalysisTaskEmcal::AddParticleContainer(const char *n)
 }
 
 //________________________________________________________________________
-AliClusterContainer* AliAnalysisTaskEmcal::AddClusterContainer(const char *n) {
-
+AliClusterContainer* AliAnalysisTaskEmcal::AddClusterContainer(const char *n) 
+{
   // Add cluster container
   // will be called in AddTask macro
 
@@ -1042,7 +1106,8 @@ AliClusterContainer* AliAnalysisTaskEmcal::AddClusterContainer(const char *n) {
 }
 
 //________________________________________________________________________
-AliParticleContainer* AliAnalysisTaskEmcal::GetParticleContainer(Int_t i) const {
+AliParticleContainer* AliAnalysisTaskEmcal::GetParticleContainer(Int_t i) const 
+{
   // Get i^th particle container
 
   if (i<0 || i>fParticleCollArray.GetEntriesFast()) return 0;
@@ -1051,7 +1116,8 @@ AliParticleContainer* AliAnalysisTaskEmcal::GetParticleContainer(Int_t i) const
 }
 
 //________________________________________________________________________
-AliClusterContainer* AliAnalysisTaskEmcal::GetClusterContainer(Int_t i) const {
+AliClusterContainer* AliAnalysisTaskEmcal::GetClusterContainer(Int_t i) const 
+{
   // Get i^th cluster container
 
   if (i<0 || i>fClusterCollArray.GetEntriesFast()) return 0;
@@ -1060,7 +1126,8 @@ AliClusterContainer* AliAnalysisTaskEmcal::GetClusterContainer(Int_t i) const {
 }
 
 //________________________________________________________________________
-AliParticleContainer* AliAnalysisTaskEmcal::GetParticleContainer(const char *name) const {
+AliParticleContainer* AliAnalysisTaskEmcal::GetParticleContainer(const char *name) const 
+{
   // Get particle container with name
 
   AliParticleContainer *cont = static_cast<AliParticleContainer*>(fParticleCollArray.FindObject(name));
@@ -1068,7 +1135,8 @@ AliParticleContainer* AliAnalysisTaskEmcal::GetParticleContainer(const char *nam
 }
 
 //________________________________________________________________________
-AliClusterContainer* AliAnalysisTaskEmcal::GetClusterContainer(const char *name) const {
+AliClusterContainer* AliAnalysisTaskEmcal::GetClusterContainer(const char *name) const 
+{
   // Get cluster container with name
 
   AliClusterContainer *cont = static_cast<AliClusterContainer*>(fClusterCollArray.FindObject(name));
@@ -1076,7 +1144,8 @@ AliClusterContainer* AliAnalysisTaskEmcal::GetClusterContainer(const char *name)
 }
 
 //________________________________________________________________________
-TClonesArray* AliAnalysisTaskEmcal::GetParticleArray(Int_t i) const {
+TClonesArray* AliAnalysisTaskEmcal::GetParticleArray(Int_t i) const 
+{
   // Get i^th TClonesArray with AliVParticle
 
   AliParticleContainer *cont = GetParticleContainer(i);
@@ -1089,7 +1158,8 @@ TClonesArray* AliAnalysisTaskEmcal::GetParticleArray(Int_t i) const {
 }
 
 //________________________________________________________________________
-TClonesArray* AliAnalysisTaskEmcal::GetClusterArray(Int_t i) const {
+TClonesArray* AliAnalysisTaskEmcal::GetClusterArray(Int_t i) const 
+{
   // Get i^th TClonesArray with AliVCluster
 
   AliClusterContainer *cont = GetClusterContainer(i);
@@ -1101,7 +1171,8 @@ TClonesArray* AliAnalysisTaskEmcal::GetClusterArray(Int_t i) const {
 }
 
 //________________________________________________________________________
-AliVParticle* AliAnalysisTaskEmcal::GetAcceptParticleFromArray(Int_t p, Int_t c) const {
+AliVParticle* AliAnalysisTaskEmcal::GetAcceptParticleFromArray(Int_t p, Int_t c) const 
+{
   // Get particle p if accepted from  container c
   // If particle not accepted return 0
 
@@ -1116,7 +1187,8 @@ AliVParticle* AliAnalysisTaskEmcal::GetAcceptParticleFromArray(Int_t p, Int_t c)
 }
 
 //________________________________________________________________________
-AliVCluster* AliAnalysisTaskEmcal::GetAcceptClusterFromArray(Int_t cl, Int_t c) const {
+AliVCluster* AliAnalysisTaskEmcal::GetAcceptClusterFromArray(Int_t cl, Int_t c) const 
+{
   // Get particle p if accepted from  container c
   // If particle not accepted return 0
 
@@ -1131,7 +1203,8 @@ AliVCluster* AliAnalysisTaskEmcal::GetAcceptClusterFromArray(Int_t cl, Int_t c)
 }
 
 //________________________________________________________________________
-Int_t AliAnalysisTaskEmcal::GetNParticles(Int_t i) const {
+Int_t AliAnalysisTaskEmcal::GetNParticles(Int_t i) const 
+{
   // Get number of entries in particle array i
 
   AliParticleContainer *cont = GetParticleContainer(i);
@@ -1143,7 +1216,8 @@ Int_t AliAnalysisTaskEmcal::GetNParticles(Int_t i) const {
 }
 
 //________________________________________________________________________
-Int_t AliAnalysisTaskEmcal::GetNClusters(Int_t i) const {
+Int_t AliAnalysisTaskEmcal::GetNClusters(Int_t i) const 
+{
   // Get number of entries in cluster array i
 
   AliClusterContainer *cont = GetClusterContainer(i);
@@ -1155,11 +1229,14 @@ Int_t AliAnalysisTaskEmcal::GetNClusters(Int_t i) const {
 }
 
 //________________________________________________________________________
-AliEmcalTriggerPatchInfo* AliAnalysisTaskEmcal::GetMainTriggerPatch() {
-  //get main trigger match; if not known yet, look for it and cache
-
-  if (fMainTriggerPatch) 
-    return fMainTriggerPatch;
+AliEmcalTriggerPatchInfo* AliAnalysisTaskEmcal::GetMainTriggerPatch(TriggerCategory trigger, Bool_t doSimpleOffline)
+{
+  // Get main trigger match
+  //
+  // For the selection of the main trigger patch, high and low threshold triggers of a given category are grouped
+  // If there are more than 1 main patch of a given trigger category (i.e. different high and low threshold patches),
+  // the highest one according to the ADC value is taken. In case doSimpleOffline is true, then only the patches from
+  // the simple offline trigger are used.
 
   if (!fTriggerPatchInfo) {
     AliError(Form("%s: fTriggerPatchInfo not available",GetName()));
@@ -1169,16 +1246,141 @@ AliEmcalTriggerPatchInfo* AliAnalysisTaskEmcal::GetMainTriggerPatch() {
   //number of patches in event
   Int_t nPatch = fTriggerPatchInfo->GetEntries();
 
-  //extract main trigger patch
-  AliEmcalTriggerPatchInfo *patch;
+  //extract main trigger patch(es)
+  AliEmcalTriggerPatchInfo *patch(NULL), *selected(NULL);
   for (Int_t iPatch = 0; iPatch < nPatch; iPatch++) {
-    
+
     patch = (AliEmcalTriggerPatchInfo*)fTriggerPatchInfo->At( iPatch );
     if (patch->IsMainTrigger()) {
-      fMainTriggerPatch = patch;
-      break;
+      if(doSimpleOffline){
+        if(patch->IsOfflineSimple()){
+          switch(trigger){
+          case kTriggerLevel0:
+            // option not yet implemented in the trigger maker
+            if(patch->IsLevel0()) selected = patch;
+            break;
+          case kTriggerLevel1Jet: 
+            if(patch->IsJetHighSimple() || patch->IsJetLowSimple()){
+              if(!selected) selected = patch;
+              else if(patch->GetADCOfflineAmp() > selected->GetADCOfflineAmp()) selected = patch;
+            }
+            break;
+          case kTriggerLevel1Gamma:
+           if(patch->IsGammaHighSimple() || patch->IsGammaLowSimple()){
+             if(!selected) selected = patch;
+             else if(patch->GetADCOfflineAmp() > selected->GetADCOfflineAmp()) selected = patch;
+           }
+           break;
+          default:   // Silence compiler warnings
+            AliError("Untreated case: Main Patch is recalculated; should be in 'else' branch");
+          };
+        }
+      } else {  // Not OfflineSimple
+        switch(trigger){
+        case kTriggerLevel0:
+            if(patch->IsLevel0()) selected = patch;
+            break;
+        case kTriggerLevel1Jet:
+            if(patch->IsJetHigh() || patch->IsJetLow()){
+              if(!selected) selected = patch;
+              else if (patch->GetADCAmp() > selected->GetADCAmp()) 
+                selected = patch;
+            }
+            break;
+        case kTriggerLevel1Gamma:
+            if(patch->IsGammaHigh() || patch->IsGammaLow()){
+              if(!selected) selected = patch;
+              else if (patch->GetADCAmp() > selected->GetADCAmp()) 
+                selected = patch;
+            }
+            break;
+         default:
+            AliError("Untreated case: Main Patch is recalculated; should be in 'else' branch");
+        };
+      }
+    }
+    else if ((trigger == kTriggerRecalcJet &&  patch->IsRecalcJet()) || 
+             (trigger == kTriggerRecalcGamma && patch->IsRecalcGamma())) {  // recalculated patches
+      if (doSimpleOffline && patch->IsOfflineSimple()) {
+        if(!selected) selected = patch;
+        else if (patch->GetADCOfflineAmp() > selected->GetADCOfflineAmp())  // this in fact should not be needed, but we have it in teh other branches as well, so keeping it for compleness
+          selected = patch;
+      }
+      else if (!doSimpleOffline && !patch->IsOfflineSimple()) {
+        if(!selected) selected = patch;
+        else if (patch->GetADCAmp() > selected->GetADCAmp()) 
+          selected = patch;
+      }
     }
   }
+  return selected;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEmcal::AddObjectToEvent(TObject *obj)
+{
+  // Add object to event
 
-  return fMainTriggerPatch;
+  if (!(InputEvent()->FindListObject(obj->GetName()))) {
+    InputEvent()->AddObject(obj);
+  } else {
+    AliFatal(Form("%s: Container with name %s already present. Aborting", GetName(), obj->GetName()));
+  }
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEmcal::GenerateFixedBinArray(Int_t n, Double_t min, Double_t max, Double_t* array) const
+{
+  Double_t binWidth = (max-min)/n;
+  array[0] = min;
+  for (Int_t i = 1; i <= n; i++) {
+    array[i] = array[i-1]+binWidth;
+  }
+}
+
+//________________________________________________________________________
+Double_t* AliAnalysisTaskEmcal::GenerateFixedBinArray(Int_t n, Double_t min, Double_t max) const
+{
+  Double_t *array = new Double_t[n+1];
+
+  GenerateFixedBinArray(n, min, max, array);
+
+  return array;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEmcal::SetRejectionReasonLabels(TAxis* axis)
+{
+  axis->SetBinLabel(1,  "NullObject");
+  axis->SetBinLabel(2,  "Pt");
+  axis->SetBinLabel(3,  "Acceptance");
+  axis->SetBinLabel(4,  "BitMap");
+  axis->SetBinLabel(5,  "Bit4");
+  axis->SetBinLabel(6,  "Bit5");
+  axis->SetBinLabel(7,  "Bit6");
+  axis->SetBinLabel(8,  "Bit7");
+  axis->SetBinLabel(9,  "MCFlag");
+  axis->SetBinLabel(10, "MCGenerator");
+  axis->SetBinLabel(11, "ChargeCut");
+  axis->SetBinLabel(12, "MinDistanceTPCSectorEdge");
+  axis->SetBinLabel(13, "MinMCLabelAccept");
+  axis->SetBinLabel(14, "IsEMCal");
+  axis->SetBinLabel(15, "Time");
+  axis->SetBinLabel(16, "Energy");
+  axis->SetBinLabel(17, "Bit16");
+  axis->SetBinLabel(18, "Bit17");
+  axis->SetBinLabel(19, "Area");
+  axis->SetBinLabel(20, "AreaEmc");
+  axis->SetBinLabel(21, "ZLeadingCh");
+  axis->SetBinLabel(22, "ZLeadingEmc");
+  axis->SetBinLabel(23, "NEF");
+  axis->SetBinLabel(24, "MinLeadPt");
+  axis->SetBinLabel(25, "MaxTrackPt");
+  axis->SetBinLabel(26, "MaxClusterPt");
+  axis->SetBinLabel(27, "Flavour");
+  axis->SetBinLabel(28, "TagStatus");
+  axis->SetBinLabel(29, "Bit28");
+  axis->SetBinLabel(30, "Bit29");
+  axis->SetBinLabel(31, "Bit30");
+  axis->SetBinLabel(32, "Bit31");
 }