]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG0/dNdEta/AlidNdEtaTask.cxx
Update
[u/mrichter/AliRoot.git] / PWG0 / dNdEta / AlidNdEtaTask.cxx
index 573ef5a4d0440c18e6180b5a0e66d7e85f8f795f..36c03509e50740915108fbd539fe9c5a2cae0386 100644 (file)
@@ -9,6 +9,7 @@
 #include <TChain.h>
 #include <TFile.h>
 #include <TH1F.h>
+#include <TH1D.h>
 #include <TH2F.h>
 #include <TH3F.h>
 #include <TParticle.h>
@@ -16,6 +17,7 @@
 #include <TNtuple.h>
 #include <TObjString.h>
 #include <TF1.h>
+#include <TGraph.h>
 
 #include <AliLog.h>
 #include <AliESDVertex.h>
 #include <AliMCEventHandler.h>
 #include <AliMCEvent.h>
 #include <AliESDInputHandler.h>
+#include <AliESDInputHandlerRP.h>
+#include <AliESDHeader.h>
 
-#include "esdTrackCuts/AliESDtrackCuts.h"
+#include "AliESDtrackCuts.h"
 #include "AliPWG0Helper.h"
 #include "AliCorrection.h"
 #include "AliCorrectionMatrix3D.h"
 #include "dNdEta/dNdEtaAnalysis.h"
+#include "AliTriggerAnalysis.h"
+#include "AliPhysicsSelection.h"
+
+//#define FULLALIROOT
+
+#ifdef FULLALIROOT
+  #include "../ITS/AliITSRecPoint.h"
+  #include "AliCDBManager.h"
+  #include "AliCDBEntry.h"
+  #include "AliGeomManager.h"
+  #include "TGeoManager.h"
+#endif
 
 ClassImp(AlidNdEtaTask)
 
 AlidNdEtaTask::AlidNdEtaTask(const char* opt) :
-  AliAnalysisTask("AlidNdEtaTask", ""),
+  AliAnalysisTaskSE("AlidNdEtaTask"),
   fESD(0),
   fOutput(0),
   fOption(opt),
-  fAnalysisMode(AliPWG0Helper::kTPC),
+  fAnalysisMode((AliPWG0Helper::AnalysisMode) (AliPWG0Helper::kTPC | AliPWG0Helper::kFieldOn)),
+  fTrigger(AliTriggerAnalysis::kMB1),
+  fFillPhi(kFALSE),
+  fDeltaPhiCut(-1),
   fReadMC(kFALSE),
   fUseMCVertex(kFALSE),
+  fOnlyPrimaries(kFALSE),
   fUseMCKine(kFALSE),
+  fSymmetrize(kFALSE),
+  fMultAxisEta1(kFALSE),
+  fDiffTreatment(AliPWG0Helper::kMCFlags),
   fEsdTrackCuts(0),
   fdNdEtaAnalysisESD(0),
   fMult(0),
@@ -53,20 +76,50 @@ AlidNdEtaTask::AlidNdEtaTask(const char* opt) :
   fEvents(0),
   fVertexResolution(0),
   fdNdEtaAnalysis(0),
+  fdNdEtaAnalysisND(0),
+  fdNdEtaAnalysisNSD(0),
+  fdNdEtaAnalysisOnePart(0),
   fdNdEtaAnalysisTr(0),
   fdNdEtaAnalysisTrVtx(0),
   fdNdEtaAnalysisTracks(0),
-  fVertex(0),
   fPartPt(0),
-  fDeltaPhi(0)
+  fVertex(0),
+  fVertexVsMult(0),
+  fPhi(0),
+  fRawPt(0),
+  fEtaPhi(0),
+  fModuleMap(0),
+  fDeltaPhi(0),
+  fDeltaTheta(0),
+  fFiredChips(0),
+  fTrackletsVsClusters(0),
+  fTrackletsVsUnassigned(0),
+  fStats(0),
+  fStats2(0),
+  fPtMin(0.15),
+  fEta(0x0),
+  fEtaMC(0x0),
+  fHistEvents(0),
+  fHistEventsMC(0),
+  fTrigEffNum(0),
+  fTrigEffDen(0),
+  fVtxEffNum(0),
+  fVtxEffDen(0),
+  fVtxTrigEffNum(0),
+  fVtxTrigEffDen(0)
 {
   //
   // Constructor. Initialization of pointers
   //
 
   // Define input and output slots here
-  DefineInput(0, TChain::Class());
-  DefineOutput(0, TList::Class());
+  //DefineInput(0, TChain::Class());
+  DefineOutput(1, TList::Class());
+  
+  fZPhi[0] = 0;
+  fZPhi[1] = 0;
+
+  AliLog::SetClassDebugLevel("AlidNdEtaTask", AliLog::kWarning);
 }
 
 AlidNdEtaTask::~AlidNdEtaTask()
@@ -78,59 +131,83 @@ AlidNdEtaTask::~AlidNdEtaTask()
   // histograms are in the output list and deleted when the output
   // list is deleted by the TSelector dtor
 
-  if (fOutput) {
+  if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
     delete fOutput;
     fOutput = 0;
   }
 }
 
+Bool_t AlidNdEtaTask::UserNotify()
+{
+  static Int_t count = 0;
+  count++;
+  Printf("Processing %d. file: %s", count, ((TTree*) GetInputData(0))->GetCurrentFile()->GetName());
+  return kTRUE;
+}
+
 //________________________________________________________________________
-void AlidNdEtaTask::ConnectInputData(Option_t *)
+void AlidNdEtaTask::ConnectInputData(Option_t *opt)
 {
   // Connect ESD
   // Called once
+  
+  AliAnalysisTaskSE::ConnectInputData(opt);
 
   Printf("AlidNdEtaTask::ConnectInputData called");
 
-  TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
-  if (!tree) {
-    Printf("ERROR: Could not read tree from input slot 0");
-  } else {
-    // Disable all branches and enable only the needed ones
-    //tree->SetBranchStatus("*", 0);
-
-    tree->SetBranchStatus("fTriggerMask", 1);
-    tree->SetBranchStatus("fSPDVertex*", 1);
-    // PrimaryVertex also needed
+  /*
+  AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 
-    if (fAnalysisMode == AliPWG0Helper::kSPD) {
-      tree->SetBranchStatus("fSPDMult*", 1);
-      //AliPWG0Helper::SetBranchStatusRecursive(tree, "AliMultiplicity", kTRUE);
-    }
+  if (!esdH) {
+    Printf("ERROR: Could not get ESDInputHandler");
+  } else {
+    fESD = esdH->GetEvent();
     
-    if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS) {
-      AliESDtrackCuts::EnableNeededBranches(tree);
-      tree->SetBranchStatus("fTracks.fLabel", 1);
-      //AliPWG0Helper::SetBranchStatusRecursive(tree, "Tracks", kTRUE);
-    }
+    TString branches("AliESDHeader Vertex ");
 
-    AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
-
-    if (!esdH) {
-      Printf("ERROR: Could not get ESDInputHandler");
-    } else
-      fESD = esdH->GetEvent();
+    if (fAnalysisMode & AliPWG0Helper::kSPD || fTrigger & AliTriggerAnalysis::kOfflineFlag)
+      branches += "AliMultiplicity ";
+      
+    if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
+      branches += "Tracks ";
+  
+    // Enable only the needed branches
+    esdH->SetActiveBranches(branches);
   }
+  */
 
   // disable info messages of AliMCEvent (per event)
   AliLog::SetClassDebugLevel("AliMCEvent", AliLog::kWarning - AliLog::kDebug + 1);
+  
+  #ifdef FULLALIROOT
+    AliCDBManager::Instance()->SetDefaultStorage("raw://");
+    //AliCDBManager::Instance()->SetDefaultStorage("MC", "Residual");
+    AliCDBManager::Instance()->SetRun(0);
+    
+    AliCDBManager* mgr = AliCDBManager::Instance();
+    AliCDBEntry* obj = mgr->Get(AliCDBPath("GRP", "Geometry", "Data"));
+    AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject());
+    
+    AliGeomManager::GetNalignable("ITS");
+    AliGeomManager::ApplyAlignObjsFromCDB("ITS");
+  #endif
 }
 
-void AlidNdEtaTask::CreateOutputObjects()
+void AlidNdEtaTask::UserCreateOutputObjects()
 {
   // create result objects and add to output list
 
   Printf("AlidNdEtaTask::CreateOutputObjects");
+  //AliLog::SetClassDebugLevel("AliPhysicsSelection", AliLog::kDebug);
+
+  if (fOnlyPrimaries)
+    Printf("WARNING: Processing only primaries (MC information used). This is for systematical checks only.");
+
+  if (fUseMCKine)
+    Printf("WARNING: Using MC kine information. This is for systematical checks only.");
+
+  if (fUseMCVertex)
+    Printf("WARNING: Replacing vertex by MC vertex. This is for systematical checks only.");
 
   fOutput = new TList;
   fOutput->SetOwner();
@@ -146,22 +223,104 @@ void AlidNdEtaTask::CreateOutputObjects()
 
   for (Int_t i=0; i<3; ++i)
   {
-    fPartEta[i] = new TH1F(Form("dndeta_check_%d", i), Form("dndeta_check_%d", i), 60, -6, 6);
+    fPartEta[i] = new TH1F(Form("dndeta_check_%d", i), Form("dndeta_check_%d", i), 60, -3, 3);
     fPartEta[i]->Sumw2();
     fOutput->Add(fPartEta[i]);
   }
 
-  fEvents = new TH1F("dndeta_check_vertex", "dndeta_check_vertex", 40, -20, 20);
+  fEvents = new TH1F("dndeta_check_vertex", "dndeta_check_vertex", 800, -40, 40);
   fOutput->Add(fEvents);
 
-  fVertexResolution = new TH1F("dndeta_vertex_resolution_z", "dndeta_vertex_resolution_z", 1000, 0, 10);
+  fVertexResolution = new TH2F("dndeta_vertex_resolution_z", ";z resolution;#Delta #phi", 1000, 0, 2, 100, 0, 0.2);
   fOutput->Add(fVertexResolution);
 
+  fPhi = new TH1F("fPhi", "fPhi;#phi in rad.;count", 720, 0, 2 * TMath::Pi());
+  fOutput->Add(fPhi);
+
+  fEtaPhi = new TH2F("fEtaPhi", "fEtaPhi;#eta;#phi in rad.;count", 80, -4, 4, 18*20, 0, 2 * TMath::Pi());
+  fOutput->Add(fEtaPhi);
+  
+  fStats = new TH1F("fStats", "fStats", 5, 0.5, 5.5);
+  fStats->GetXaxis()->SetBinLabel(1, "vertexer 3d");
+  fStats->GetXaxis()->SetBinLabel(2, "vertexer z");
+  fStats->GetXaxis()->SetBinLabel(3, "trigger");
+  fStats->GetXaxis()->SetBinLabel(4, "physics events");
+  fStats->GetXaxis()->SetBinLabel(5, "physics events after veto");
+  fOutput->Add(fStats);
+  
+  fStats2 = new TH2F("fStats2", "fStats2", 7, -0.5, 6.5, 12, -0.5, 11.5);
+  
+  fStats2->GetXaxis()->SetBinLabel(1, "No trigger");
+  fStats2->GetXaxis()->SetBinLabel(2, "No Vertex");
+  fStats2->GetXaxis()->SetBinLabel(3, "Vtx res. too bad");
+  fStats2->GetXaxis()->SetBinLabel(4, "|z-vtx| > 15");
+  fStats2->GetXaxis()->SetBinLabel(5, "|z-vtx| > 10");
+  fStats2->GetXaxis()->SetBinLabel(6, "0 tracklets");
+  fStats2->GetXaxis()->SetBinLabel(7, "Selected");
+  
+  fStats2->GetYaxis()->SetBinLabel(1, "n/a");
+  fStats2->GetYaxis()->SetBinLabel(2, "empty");
+  fStats2->GetYaxis()->SetBinLabel(3, "BBA");
+  fStats2->GetYaxis()->SetBinLabel(4, "BBC");
+  fStats2->GetYaxis()->SetBinLabel(5, "BBA BBC");
+  fStats2->GetYaxis()->SetBinLabel(6, "BGA");
+  fStats2->GetYaxis()->SetBinLabel(7, "BGC");
+  fStats2->GetYaxis()->SetBinLabel(8, "BGA BGC");
+  fStats2->GetYaxis()->SetBinLabel(9, "BBA BGC");
+  fStats2->GetYaxis()->SetBinLabel(10, "BGA BBC");
+  fStats2->GetYaxis()->SetBinLabel(11, "GFO");
+  fStats2->GetYaxis()->SetBinLabel(12, "NO GFO");
+  fOutput->Add(fStats2);
+
+  fTrackletsVsClusters = new TH2F("fTrackletsVsClusters", ";tracklets;clusters in ITS", 50, -0.5, 49.5, 1000, -0.5, 999.5);
+  fOutput->Add(fTrackletsVsClusters);
+  
+  if (fAnalysisMode & AliPWG0Helper::kSPD)
+  {
+    fDeltaPhi = new TH1F("fDeltaPhi", "fDeltaPhi;#Delta #phi;Entries", 500, -0.16, 0.16);
+    fOutput->Add(fDeltaPhi);
+    fDeltaTheta = new TH1F("fDeltaTheta", "fDeltaTheta;#Delta #theta;Entries", 500, -0.05, 0.05);
+    fOutput->Add(fDeltaTheta);
+    fFiredChips = new TH2F("fFiredChips", "fFiredChips;Chips L1 + L2;tracklets", 1201, -0.5, 1201, 50, -0.5, 49.5);
+    fOutput->Add(fFiredChips);
+    fTrackletsVsUnassigned = new TH2F("fTrackletsVsUnassigned", ";tracklets;unassigned clusters in L0", 50, -0.5, 49.5, 200, -0.5, 199.5);
+    fOutput->Add(fTrackletsVsUnassigned);
+    for (Int_t i=0; i<2; i++)
+    {
+      fZPhi[i] = new TH2F(Form("fZPhi_%d", i), Form("fZPhi Layer %d;z (cm);#phi (rad.)", i), 200, -15, 15, 200, 0, TMath::Pi() * 2);
+      fOutput->Add(fZPhi[i]);
+    }
+    
+    fModuleMap = new TH1F("fModuleMap", "fModuleMap;module number;cluster count", 240, -0.5, 239.5);
+    fOutput->Add(fModuleMap);
+  }
+
+  if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
+  {
+    fRawPt =  new TH1F("fRawPt", "raw pt;p_{T};Count", 2000, 0, 100);
+    fOutput->Add(fRawPt);
+  }
+
+  fVertex = new TH3F("vertex_check", "vertex_check;x;y;z", 100, -1, 1, 100, -1, 1, 100, -30, 30);
+  fOutput->Add(fVertex);
+  
+  fVertexVsMult = new TH3F("fVertexVsMult", "fVertexVsMult;x;y;multiplicity", 100, -1, 1, 100, -1, 1, 100, -0.5, 99.5);
+  fOutput->Add(fVertexVsMult);
+
   if (fReadMC)
   {
     fdNdEtaAnalysis = new dNdEtaAnalysis("dndeta", "dndeta", fAnalysisMode);
     fOutput->Add(fdNdEtaAnalysis);
 
+    fdNdEtaAnalysisND = new dNdEtaAnalysis("dndetaND", "dndetaND", fAnalysisMode);
+    fOutput->Add(fdNdEtaAnalysisND);
+
+    fdNdEtaAnalysisNSD = new dNdEtaAnalysis("dndetaNSD", "dndetaNSD", fAnalysisMode);
+    fOutput->Add(fdNdEtaAnalysisNSD);
+
+    fdNdEtaAnalysisOnePart = new dNdEtaAnalysis("dndetaOnePart", "dndetaOnePart", fAnalysisMode);
+    fOutput->Add(fdNdEtaAnalysisOnePart);
+
     fdNdEtaAnalysisTr = new dNdEtaAnalysis("dndetaTr", "dndetaTr", fAnalysisMode);
     fOutput->Add(fdNdEtaAnalysisTr);
 
@@ -171,226 +330,738 @@ void AlidNdEtaTask::CreateOutputObjects()
     fdNdEtaAnalysisTracks = new dNdEtaAnalysis("dndetaTracks", "dndetaTracks", fAnalysisMode);
     fOutput->Add(fdNdEtaAnalysisTracks);
 
-    fVertex = new TH3F("vertex_check", "vertex_check", 50, -50, 50, 50, -50, 50, 50, -50, 50);
-    fOutput->Add(fVertex);
-
     fPartPt =  new TH1F("dndeta_check_pt", "dndeta_check_pt", 1000, 0, 10);
     fPartPt->Sumw2();
     fOutput->Add(fPartPt);
   }
 
-  if (fAnalysisMode == AliPWG0Helper::kSPD)
-    fDeltaPhi = new TH1F("fDeltaPhi", "fDeltaPhi;#Delta #phi;Entries", 1000, -3.14, 3.14);
+  if (fEsdTrackCuts)
+  {
+    fEsdTrackCuts->SetName("fEsdTrackCuts");
+    fOutput->Add(fEsdTrackCuts);
+  }
+
+  fEta = new TH1D("fEta", "Eta;#eta;count", 80, -4, 4);
+  fOutput->Add(fEta);
+  
+  fEtaMC = new TH1D("fEtaMC", "Eta, MC;#eta;count", 80, -4, 4);
+  fOutput->Add(fEtaMC);
+
+  fHistEvents = new TH1D("fHistEvents", "N. of Events;accepted;count", 2, 0, 2);
+  fOutput->Add(fHistEvents);
+  
+  fHistEventsMC = new TH1D("fHistEventsMC", "N. of MC Events;accepted;count", 2, 0, 2);
+  fOutput->Add(fHistEventsMC);
+
+  fTrigEffNum = new TH1D("fTrigEffNum", "N. of triggered events", 100,0,100); 
+  fOutput->Add(fTrigEffNum);
+  fTrigEffDen = new TH1D("fTrigEffDen", "N. of true events", 100,0,100); 
+  fOutput->Add(fTrigEffDen);
+  fVtxEffNum = new TH1D("fVtxEffNum", "N. of events with vtx", 100,0,100); 
+  fOutput->Add(fVtxEffNum);
+  fVtxEffDen = new TH1D("fVtxEffDen", "N. of true events", 100,0,100); 
+  fOutput->Add(fVtxEffDen);
+  fVtxTrigEffNum = new TH1D("fVtxTrigEffNum", "N. of triggered events with vtx", 100,0,100); 
+  fOutput->Add(fVtxTrigEffNum);
+  fVtxTrigEffDen = new TH1D("fVtxTrigEffDen", "N. of triggered true events", 100,0,100); 
+  fOutput->Add(fVtxTrigEffDen);
+  
+  //AliLog::SetClassDebugLevel("AliPhysicsSelection", AliLog::kDebug);
 }
 
-void AlidNdEtaTask::Exec(Option_t*)
+Bool_t AlidNdEtaTask::IsEventInBinZero()
 {
-  // process the event
-
-  // Check prerequisites
-  if (!fESD)
+  // checks if the event goes to the 0 bin
+  
+  Bool_t isZeroBin = kTRUE;
+  fESD = (AliESDEvent*) fInputEvent;
+  
+  AliInputEventHandler* inputHandler = static_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
+  if (!inputHandler)
   {
-    AliDebug(AliLog::kError, "ESD branch not available");
-    return;
+    Printf("ERROR: Could not receive input handler");
+    return kFALSE;
+  }
+    
+  static AliTriggerAnalysis* triggerAnalysis = 0;
+  if (!triggerAnalysis)
+  {
+    AliPhysicsSelection* physicsSelection = dynamic_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
+    if (physicsSelection)
+      triggerAnalysis = physicsSelection->GetTriggerAnalysis();
   }
-
-  // trigger definition
-  Bool_t eventTriggered = AliPWG0Helper::IsEventTriggered(fESD->GetTriggerMask(), AliPWG0Helper::kMB2);
-
-  // get the ESD vertex
-  const AliESDVertex* vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
   
-  Double_t vtx[3];
-  if (vtxESD) {
-    vtxESD->GetXYZ(vtx);
+  if (!triggerAnalysis)
+  {
+    Printf("ERROR: Could not receive trigger analysis object");
+    return kFALSE;
   }
-  else
-    Printf("No vertex");
   
-  // fill z vertex resolution
-  if (vtxESD)
-    fVertexResolution->Fill(vtxESD->GetZRes());
+  if (!triggerAnalysis->IsTriggerFired(fESD, fTrigger))
+    return kFALSE;
+  
+  // 0 bin check - from Michele
+  
+  const AliMultiplicity* mult = fESD->GetMultiplicity();
+  if (!mult){
+    Printf("AlidNdEtaTask::IsBinZero: Can't get multiplicity object from ESDs");
+    return kFALSE;
+  }
+  Int_t ntracklet = mult->GetNumberOfTracklets();
+  const AliESDVertex * vtxESD = fESD->GetPrimaryVertexSPD();
+  if(vtxESD) {
+         // If there is a vertex from vertexer z with delta phi > 0.02 we
+         // don't consider it rec (we keep the event in bin0). If quality
+         // is good eneough we check the number of tracklets
+         // if the vertex is more than 15 cm away, this is autamatically bin0
+         if( TMath::Abs(vtxESD->GetZ()) <= 15 ) {
+                 if (vtxESD->IsFromVertexerZ()) {
+                         if (vtxESD->GetDispersion()<=0.02 ) {
+                                 if(ntracklet>0) isZeroBin = kFALSE;
+                         }
+                 } else if(ntracklet>0) isZeroBin = kFALSE; // if the event is not from Vz we chek the n of tracklets
+         } 
+  }
+  return isZeroBin;
+}
 
-  // post the data already here
-  PostData(0, fOutput);
+void AlidNdEtaTask::UserExec(Option_t*)
+{
+  // process the event
 
-  // needed for syst. studies
-  AliStack* stack = 0;
+  fESD = (AliESDEvent*) fInputEvent;
   
-  if (fUseMCVertex || fUseMCKine) {
-    if (!fReadMC) {
-      Printf("ERROR: fUseMCVertex or fUseMCKine set without fReadMC set!");
-      return;
-    }
+  // these variables are also used in the MC section below; however, if ESD is off they stay with their default values
+  Bool_t eventTriggered = kFALSE;
+  const AliESDVertex* vtxESD = 0;
 
-    AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
-    if (!eventHandler) {
-      Printf("ERROR: Could not retrieve MC event handler");
-      return;
-    }
+  // post the data already here
+  PostData(1, fOutput);
 
-    AliMCEvent* mcEvent = eventHandler->MCEvent();
-    if (!mcEvent) {
-      Printf("ERROR: Could not retrieve MC event");
+  // ESD analysis
+  if (fESD)
+  {
+    AliESDHeader* esdHeader = fESD->GetHeader();
+    if (!esdHeader)
+    {
+      Printf("ERROR: esdHeader could not be retrieved");
       return;
     }
-
-    AliHeader* header = mcEvent->Header();
-    if (!header)
+    
+    AliInputEventHandler* inputHandler = (AliInputEventHandler*) AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
+    if (!inputHandler)
     {
-      AliDebug(AliLog::kError, "Header not available");
+      Printf("ERROR: Could not receive input handler");
       return;
     }
-
-    if (fUseMCVertex)
+    
+    // TODO use flags here!
+    eventTriggered = inputHandler->IsEventSelected();
+        
+    static AliTriggerAnalysis* triggerAnalysis = 0;
+    if (!triggerAnalysis)
     {
-      Printf("WARNING: Replacing vertex by MC vertex. This is for systematical checks only.");
-      // get the MC vertex
-      AliGenEventHeader* genHeader = header->GenEventHeader();
-      TArrayF vtxMC(3);
-      genHeader->PrimaryVertex(vtxMC);
-
-      vtx[2] = vtxMC[2];
+      AliPhysicsSelection* physicsSelection = dynamic_cast<AliPhysicsSelection*> (inputHandler->GetEventSelection());
+      if (physicsSelection)
+        triggerAnalysis = physicsSelection->GetTriggerAnalysis();
     }
-
-    if (fUseMCKine)
+      
+    if (eventTriggered)
+      eventTriggered = triggerAnalysis->IsTriggerFired(fESD, fTrigger);
+    
+    AliTriggerAnalysis::V0Decision v0A = triggerAnalysis->V0Trigger(fESD, AliTriggerAnalysis::kASide, kFALSE);
+    AliTriggerAnalysis::V0Decision v0C = triggerAnalysis->V0Trigger(fESD, AliTriggerAnalysis::kCSide, kFALSE);
+    Bool_t fastORHW = (triggerAnalysis->SPDFiredChips(fESD, 1) > 0);
+    
+    Int_t vZero = 0;
+    if (v0A != AliTriggerAnalysis::kV0Invalid && v0C != AliTriggerAnalysis::kV0Invalid)
     {
-
-      stack = mcEvent->Stack();
-      if (!stack)
-      {
-        AliDebug(AliLog::kError, "Stack not available");
-        return;
-      }
+      if (v0A == AliTriggerAnalysis::kV0Empty && v0C == AliTriggerAnalysis::kV0Empty) vZero = 1;
+      if (v0A == AliTriggerAnalysis::kV0BB    && v0C == AliTriggerAnalysis::kV0Empty) vZero = 2;
+      if (v0A == AliTriggerAnalysis::kV0Empty && v0C == AliTriggerAnalysis::kV0BB)    vZero = 3;
+      if (v0A == AliTriggerAnalysis::kV0BB    && v0C == AliTriggerAnalysis::kV0BB)    vZero = 4;
+      if (v0A == AliTriggerAnalysis::kV0BG    && v0C == AliTriggerAnalysis::kV0Empty) vZero = 5;
+      if (v0A == AliTriggerAnalysis::kV0Empty && v0C == AliTriggerAnalysis::kV0BG)    vZero = 6;
+      if (v0A == AliTriggerAnalysis::kV0BG    && v0C == AliTriggerAnalysis::kV0BG)    vZero = 7;
+      if (v0A == AliTriggerAnalysis::kV0BB    && v0C == AliTriggerAnalysis::kV0BG)    vZero = 8;
+      if (v0A == AliTriggerAnalysis::kV0BG    && v0C == AliTriggerAnalysis::kV0BB)    vZero = 9;
     }
-  }
 
-  // create list of (label, eta, pt) tuples
-  Int_t inputCount = 0;
-  Int_t* labelArr = 0;
-  Float_t* etaArr = 0;
-  Float_t* ptArr = 0;
-  if (fAnalysisMode == AliPWG0Helper::kSPD)
-  {
-    // get tracklets
-    const AliMultiplicity* mult = fESD->GetMultiplicity();
-    if (!mult)
+    if (vZero == 0)
+      Printf("VZERO: %d %d %d %d", vZero, eventTriggered, v0A, v0C);
+      
+    Bool_t filled = kFALSE;
+      
+    if (!eventTriggered)
     {
-      AliDebug(AliLog::kError, "AliMultiplicity not available");
-      return;
+      fStats2->Fill(0.0, vZero);
+      fStats2->Fill(0.0, (fastORHW) ? 10 : 11);
+      filled = kTRUE;
     }
+    
+    if (eventTriggered)
+      fStats->Fill(3);
+      
+    /*
+    // ITS cluster tree
+    AliESDInputHandlerRP* handlerRP = dynamic_cast<AliESDInputHandlerRP*> (AliAnalysisManager::GetAnalys
+isManager()->GetInputEventHandler());
+    if (handlerRP)
+    {
+      TTree* itsClusterTree = handlerRP->GetTreeR("ITS");
+      if (!itsClusterTree)
+        return;
 
-    labelArr = new Int_t[mult->GetNumberOfTracklets()];
-    etaArr = new Float_t[mult->GetNumberOfTracklets()];
-    ptArr = new Float_t[mult->GetNumberOfTracklets()];
+      TClonesArray* itsClusters = new TClonesArray("AliITSRecPoint");
+      TBranch* itsClusterBranch=itsClusterTree->GetBranch("ITSRecPoints");
 
-    if (fUseMCKine && stack)
-      Printf("Processing only primaries (MC information used). This is for systematical checks only.");
+      itsClusterBranch->SetAddress(&itsClusters);
 
-    // get multiplicity from ITS tracklets
-    for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
-    {
-      //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
+      Int_t nItsSubs = (Int_t)itsClusterTree->GetEntries();
+
+      Int_t totalClusters = 0;
 
-      if (fUseMCKine && stack)
-        if (mult->GetLabel(i, 0) < 0 || mult->GetLabel(i, 0) != mult->GetLabel(i, 1) || !stack->IsPhysicalPrimary(mult->GetLabel(i, 0)))
+      // loop over the its subdetectors
+      for (Int_t iIts=0; iIts < nItsSubs; iIts++) {
+
+        if (!itsClusterTree->GetEvent(iIts))
           continue;
-      
-      Float_t deltaPhi = mult->GetDeltaPhi(i);
-      // prevent values to be shifted by 2 Pi()
-      if (deltaPhi < -TMath::Pi())
-        deltaPhi += TMath::Pi() * 2;
-      if (deltaPhi > TMath::Pi())
-        deltaPhi -= TMath::Pi() * 2;
 
-      if (TMath::Abs(deltaPhi) > 1)
-        printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), deltaPhi);
+        Int_t nClusters = itsClusters->GetEntriesFast();
+        totalClusters += nClusters;
 
-      fDeltaPhi->Fill(deltaPhi);
+        #ifdef FULLALIROOT
+          if (fAnalysisMode & AliPWG0Helper::kSPD)
+          {
+            // loop over clusters
+            while (nClusters--) {
+              AliITSRecPoint* cluster = (AliITSRecPoint*) itsClusters->UncheckedAt(nClusters);
 
-      etaArr[inputCount] = mult->GetEta(i);
-      labelArr[inputCount] = mult->GetLabel(i, 0);
-      ptArr[inputCount] = 0; // no pt for tracklets
-      ++inputCount;
-    }
+              Int_t layer = cluster->GetLayer();
 
-    Printf("Accepted %d tracklets", inputCount);
-  }
-  else if (fAnalysisMode == AliPWG0Helper::kTPC || fAnalysisMode == AliPWG0Helper::kTPCITS)
-  {
-    if (!fEsdTrackCuts)
-    {
-      AliDebug(AliLog::kError, "fESDTrackCuts not available");
-      return;
-    }
+              if (layer > 1)
+                continue;
 
-    // get multiplicity from ESD tracks
-    TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD);
-    Int_t nGoodTracks = list->GetEntries();
+              Float_t xyz[3] = {0., 0., 0.};
+              cluster->GetGlobalXYZ(xyz);
 
-    labelArr = new Int_t[nGoodTracks];
-    etaArr = new Float_t[nGoodTracks];
-    ptArr = new Float_t[nGoodTracks];
+              Float_t phi = TMath::Pi() + TMath::ATan2(-xyz[1], -xyz[0]);
+              Float_t z = xyz[2];
 
-    // loop over esd tracks
-    for (Int_t i=0; i<nGoodTracks; i++)
+              fZPhi[layer]->Fill(z, phi);
+              fModuleMap->Fill(layer * 80 + cluster->GetDetectorIndex());
+            }
+          }
+        #endif
+      }
+    }
+    */
+      
+    vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
+    if (!vtxESD)
     {
-      AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
-      if (!esdTrack)
+      if (!filled)
       {
-        AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
-        continue;
+        fStats2->Fill(1, vZero);
+        fStats2->Fill(1, (fastORHW) ? 10 : 11);
+        filled = kTRUE;
+      }
+    }
+    else
+    {
+      if (!AliPWG0Helper::TestVertex(vtxESD, fAnalysisMode))
+      {
+        if (!filled)
+        {
+          fStats2->Fill(2, vZero);
+          fStats2->Fill(2, (fastORHW) ? 10 : 11);
+          filled = kTRUE;
+        }
+      }
+      
+      Double_t vtx[3];
+      vtxESD->GetXYZ(vtx);
+      
+      // try to compare spd vertex and vertexer from tracks
+      // remove vertices outside +- 15 cm
+      if (TMath::Abs(vtx[2]) > 15)
+      {
+        if (!filled)
+        {
+          fStats2->Fill(3, vZero);
+          fStats2->Fill(3, (fastORHW) ? 10 : 11);
+          filled = kTRUE;
+        }
       }
+      
+      if (TMath::Abs(vtx[2]) > 10)
+      {
+        if (!filled)
+        {
+          fStats2->Fill(4, vZero);
+          fStats2->Fill(4, (fastORHW) ? 10 : 11);
+          filled = kTRUE;
+        }
+      }
+        
+      const AliMultiplicity* mult = fESD->GetMultiplicity();
+      if (!mult)
+      {
+       Printf("Returning, no Multiplicity found");
+        return;
+      }
+      
+      if (mult->GetNumberOfTracklets() == 0)
+      {
+        if (!filled)
+        {
+          fStats2->Fill(5, vZero);
+          fStats2->Fill(5, (fastORHW) ? 10 : 11);
+          filled = kTRUE;
+        }
+      }
+    }
+    
+    if (!filled)
+    {
+      fStats2->Fill(6, vZero);
+      fStats2->Fill(6, (fastORHW) ? 10 : 11);
+      //Printf("File: %s, IEV: %d, TRG: ---, Orbit: 0x%x, Period: %d, BC: %d", ((TTree*) GetInputData(0))->GetCurrentFile()->GetName(), fESD->GetEventNumberInFile(), fESD->GetOrbitNumber(),fESD->GetPeriodNumber(),fESD->GetBunchCrossNumber());
+    }
+      
+    if (eventTriggered)
+      fStats->Fill(3);
+    
+    fStats->Fill(5);
+    
+    // get the ESD vertex
+    vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
+    
+    Double_t vtx[3];
 
-      etaArr[inputCount] = esdTrack->Eta();
-      labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
-      ptArr[inputCount] = esdTrack->Pt();
-      ++inputCount;
+    // fill z vertex resolution
+    if (vtxESD)
+    {
+      if (strcmp(vtxESD->GetTitle(), "vertexer: Z") == 0)
+        fVertexResolution->Fill(vtxESD->GetZRes(), vtxESD->GetDispersion());
+      else
+        fVertexResolution->Fill(vtxESD->GetZRes(), 0);
+      
+      //if (strcmp(vtxESD->GetTitle(), "vertexer: 3D") == 0)
+      {
+        fVertex->Fill(vtxESD->GetXv(), vtxESD->GetYv(), vtxESD->GetZv());
+      }
+      
+      if (AliPWG0Helper::TestVertex(vtxESD, fAnalysisMode))
+      {
+          vtxESD->GetXYZ(vtx);
+
+          // vertex stats
+          if (strcmp(vtxESD->GetTitle(), "vertexer: 3D") == 0)
+            fStats->Fill(1);
+          
+          if (strcmp(vtxESD->GetTitle(), "vertexer: Z") == 0)
+            fStats->Fill(2);
+          
+          // remove vertices outside +- 15 cm
+          if (TMath::Abs(vtx[2]) > 15)
+            vtxESD = 0;
+      }
+      else
+        vtxESD = 0;
     }
 
-    Printf("Accepted %d tracks", nGoodTracks);
+    // needed for syst. studies
+    AliStack* stack = 0;
+    TArrayF vtxMC(3);
 
-    delete list;
-  }
-  else
-    return;
+    if (fUseMCVertex || fUseMCKine || fOnlyPrimaries || fReadMC) {
+      if (!fReadMC) {
+        Printf("ERROR: fUseMCVertex or fUseMCKine or fOnlyPrimaries set without fReadMC set!");
+        return;
+      }
 
-  // Processing of ESD information (always)
-  if (eventTriggered)
-  {
-    // control hist
-    fMult->Fill(inputCount);
-    fdNdEtaAnalysisESD->FillTriggeredEvent(inputCount);
+      AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
+      if (!eventHandler) {
+        Printf("ERROR: Could not retrieve MC event handler");
+        return;
+      }
 
-    if (vtxESD)
-    {
-      // control hist
-      fMultVtx->Fill(inputCount);
+      AliMCEvent* mcEvent = eventHandler->MCEvent();
+      if (!mcEvent) {
+        Printf("ERROR: Could not retrieve MC event");
+        return;
+      }
+
+      AliHeader* header = mcEvent->Header();
+      if (!header)
+      {
+        AliDebug(AliLog::kError, "Header not available");
+        return;
+      }
 
-      for (Int_t i=0; i<inputCount; ++i)
+      // get the MC vertex
+      AliGenEventHeader* genHeader = header->GenEventHeader();
+      if (!genHeader)
       {
-        Float_t eta = etaArr[i];
-        Float_t pt  = ptArr[i];
+        AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
+        return;
+      }
+      genHeader->PrimaryVertex(vtxMC);
 
-        fdNdEtaAnalysisESD->FillTrack(vtx[2], eta, pt);
+      if (fUseMCVertex)
+        vtx[2] = vtxMC[2];
 
-        if (TMath::Abs(vtx[2]) < 20)
+      stack = mcEvent->Stack();
+      if (!stack)
+      {
+        AliDebug(AliLog::kError, "Stack not available");
+        return;
+      }
+    }
+    
+    // create list of (label, eta, pt) tuples
+    Int_t inputCount = 0;
+    Int_t* labelArr = 0;
+    Float_t* etaArr = 0;
+    Float_t* thirdDimArr = 0;
+    if (fAnalysisMode & AliPWG0Helper::kSPD)
+    {
+      if (vtxESD)
+      {
+        // get tracklets
+        const AliMultiplicity* mult = fESD->GetMultiplicity();
+        if (!mult)
         {
-          fPartEta[0]->Fill(eta);
-
-          if (vtx[2] < 0)
-            fPartEta[1]->Fill(eta);
-          else
-            fPartEta[2]->Fill(eta);
+          AliDebug(AliLog::kError, "AliMultiplicity not available");
+          return;
+        }
+  
+        Int_t arrayLength = mult->GetNumberOfTracklets();
+        if (fAnalysisMode & AliPWG0Helper::kSPDOnlyL0)
+          arrayLength += mult->GetNumberOfSingleClusters();
+          
+        labelArr = new Int_t[arrayLength];
+        etaArr = new Float_t[arrayLength];
+        thirdDimArr = new Float_t[arrayLength];
+  
+        // get multiplicity from SPD tracklets
+        for (Int_t i=0; i<mult->GetNumberOfTracklets(); ++i)
+        {
+          //printf("%d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), mult->GetDeltaPhi(i));
+  
+          if (fOnlyPrimaries)
+            if (mult->GetLabel(i, 0) < 0 || mult->GetLabel(i, 0) != mult->GetLabel(i, 1) || !stack->IsPhysicalPrimary(mult->GetLabel(i, 0)))
+              continue;
+  
+          Float_t deltaPhi = mult->GetDeltaPhi(i);
+          // prevent values to be shifted by 2 Pi()
+          if (deltaPhi < -TMath::Pi())
+            deltaPhi += TMath::Pi() * 2;
+          if (deltaPhi > TMath::Pi())
+            deltaPhi -= TMath::Pi() * 2;
+  
+          if (TMath::Abs(deltaPhi) > 1)
+            printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, mult->GetTheta(i), mult->GetPhi(i), deltaPhi);
+  
+          Int_t label = mult->GetLabel(i, 0);
+          Float_t eta = mult->GetEta(i);
+          
+          // control histograms
+          Float_t phi = mult->GetPhi(i);
+          if (phi < 0)
+            phi += TMath::Pi() * 2;
+            
+          // TEST exclude probably inefficient phi region
+          //if (phi > 5.70 || phi < 0.06)
+          //  continue;
+            
+          fPhi->Fill(phi);
+          
+          if (vtxESD && TMath::Abs(vtx[2]) < 10)
+          {
+            fEtaPhi->Fill(eta, phi);
+            fDeltaPhi->Fill(deltaPhi);
+            fDeltaTheta->Fill(mult->GetDeltaTheta(i));
+          }
+          
+          if (fDeltaPhiCut > 0 && (TMath::Abs(deltaPhi) > fDeltaPhiCut || TMath::Abs(mult->GetDeltaTheta(i)) > fDeltaPhiCut / 0.08 * 0.025))
+            continue;
+  
+          if (fUseMCKine)
+          {
+            if (label > 0)
+            {
+              TParticle* particle = stack->Particle(label);
+              eta = particle->Eta();
+              phi = particle->Phi();
+            }
+            else
+              Printf("WARNING: fUseMCKine set without fOnlyPrimaries and no label found");
+          }
+          
+          if (fSymmetrize)
+            eta = TMath::Abs(eta);
+  
+          etaArr[inputCount] = eta;
+          labelArr[inputCount] = label;
+          thirdDimArr[inputCount] = phi;
+          ++inputCount;
+        }
+        
+        if (fAnalysisMode & AliPWG0Helper::kSPDOnlyL0)
+        {
+          // get additional clusters from L0 
+          for (Int_t i=0; i<mult->GetNumberOfSingleClusters(); ++i)
+          {
+            etaArr[inputCount] = -TMath::Log(TMath::Tan(mult->GetThetaSingle(i)/2.));
+            labelArr[inputCount] = -1;
+            thirdDimArr[inputCount] = mult->GetPhiSingle(i);
+            
+            ++inputCount;
+          }
+        }
+        
+        if (!fFillPhi)
+        {
+          // fill multiplicity in third axis
+          for (Int_t i=0; i<inputCount; ++i)
+            thirdDimArr[i] = inputCount;
         }
+  
+        Int_t firedChips = mult->GetNumberOfFiredChips(0) + mult->GetNumberOfFiredChips(1);
+        fFiredChips->Fill(firedChips, inputCount);
+        Printf("Accepted %d tracklets (%d fired chips)", inputCount, firedChips);
+        
+        fTrackletsVsUnassigned->Fill(inputCount, mult->GetNumberOfSingleClusters());
+      }
+    }
+    else if (fAnalysisMode & AliPWG0Helper::kTPC || fAnalysisMode & AliPWG0Helper::kTPCITS)
+    {
+      if (!fEsdTrackCuts)
+      {
+        AliDebug(AliLog::kError, "fESDTrackCuts not available");
+        return;
       }
 
-      // for event count per vertex
-      fdNdEtaAnalysisESD->FillEvent(vtx[2], inputCount);
+      Bool_t foundInEta10 = kFALSE;
+      
+      if (vtxESD)
+      {
+        // get multiplicity from ESD tracks
+        TObjArray* list = fEsdTrackCuts->GetAcceptedTracks(fESD, fAnalysisMode & AliPWG0Helper::kTPC);
+        Int_t nGoodTracks = list->GetEntries();
+        Printf("Accepted %d tracks out of %d total ESD tracks", nGoodTracks, fESD->GetNumberOfTracks());
+  
+        labelArr = new Int_t[nGoodTracks];
+        etaArr = new Float_t[nGoodTracks];
+        thirdDimArr = new Float_t[nGoodTracks];
+  
+        // loop over esd tracks
+        for (Int_t i=0; i<nGoodTracks; i++)
+        {
+          AliESDtrack* esdTrack = dynamic_cast<AliESDtrack*> (list->At(i));
+          if (!esdTrack)
+          {
+            AliDebug(AliLog::kError, Form("ERROR: Could not retrieve track %d.", i));
+            continue;
+          }
+          
+          Float_t phi = esdTrack->Phi();
+          if (phi < 0)
+            phi += TMath::Pi() * 2;
+  
+          Float_t eta = esdTrack->Eta();
+          Int_t label = TMath::Abs(esdTrack->GetLabel());
+          Float_t pT  = esdTrack->Pt();
+          
+          // force pT to fixed value without B field
+          if (!(fAnalysisMode & AliPWG0Helper::kFieldOn))
+            pT = 1;
+  
+          fPhi->Fill(phi);
+          fEtaPhi->Fill(eta, phi);
+          if (eventTriggered && vtxESD)
+            fRawPt->Fill(pT);
+  
+          if (esdTrack->Pt() < fPtMin) // even if the pt cut is already in the esdTrackCuts....
+            continue;
+          
+          if (fOnlyPrimaries)
+          {
+          if (label == 0)
+            continue;
+          
+          if (stack->IsPhysicalPrimary(label) == kFALSE)
+            continue;
+          }
+
+         // 2 types of INEL>0 trigger - choose one
+
+         // 1. HL 
+          // INEL>0 trigger
+          // if (TMath::Abs(esdTrack->Eta()) < 1 && esdTrack->Pt() > 0.15)
+          // foundInEta10 = kTRUE;
+
+          // 2. MB working group
+         if (TMath::Abs(esdTrack->Eta()) < 0.8 && esdTrack->Pt() > fPtMin){  // this should be in the trigger selection as well, so one particle should always be found
+           foundInEta10 = kTRUE;
+         }
+  
+          if (fUseMCKine)
+          {
+           if (label > 0)
+            {
+              TParticle* particle = stack->Particle(label);
+              eta = particle->Eta();
+              pT = particle->Pt();
+             // check when using INEL>0, MB Working Group definition
+             if (TMath::Abs(eta) >=0.8 || pT <= fPtMin){
+               AliDebug(2,Form("WARNING ************* USING KINE: eta = %f, pt = %f",eta,pT));
+             }
+            }
+            else
+              Printf("WARNING: fUseMCKine set without fOnlyPrimaries and no label found");
+         }
+  
+          if (fSymmetrize)
+            eta = TMath::Abs(eta);
+          etaArr[inputCount] = eta;
+          labelArr[inputCount] = TMath::Abs(esdTrack->GetLabel());
+          thirdDimArr[inputCount] = pT;
+          ++inputCount;
+        }
+        
+        //if (inputCount > 30)
+        //  Printf("Event with %d accepted TPC tracks. Period number: %d Orbit number: %x Bunch crossing number: %d", inputCount, fESD->GetPeriodNumber(), fESD->GetOrbitNumber(), fESD->GetBunchCrossNumber());
+        
+        // TODO restrict inputCount used as measure for the multiplicity to |eta| < 1
+  
+        delete list;
+      }
+      
+      if (!foundInEta10){
+        eventTriggered = kFALSE;
+       fHistEvents->Fill(0.1);
+       AliDebug(3,"Event rejected");
+      }
+      else{
+       if (eventTriggered){
+         fHistEvents->Fill(1.1);
+         AliDebug(3,Form("Event Accepted, with inputcount = %d",inputCount));
+       }
+       else{
+         AliDebug(3,"Event has foundInEta10 but was not triggered");
+       }
+      }
+    }
+    else
+      return;
 
+    // Processing of ESD information (always)
+    if (eventTriggered)
+    {
       // control hist
-      fEvents->Fill(vtx[2]);
+      fMult->Fill(inputCount);
+      fdNdEtaAnalysisESD->FillTriggeredEvent(inputCount);
+
+      if (vtxESD)
+      {
+        // control hist
+        
+        if (strcmp(vtxESD->GetTitle(), "vertexer: 3D") == 0)
+          fVertexVsMult->Fill(vtxESD->GetXv(), vtxESD->GetYv(), inputCount);
+      
+        Int_t part05 = 0;
+        Int_t part10 = 0;
+        for (Int_t i=0; i<inputCount; ++i)
+        {
+          Float_t eta = etaArr[i];
+          Float_t thirdDim = thirdDimArr[i];
+         fEta->Fill(eta);
+
+          fdNdEtaAnalysisESD->FillTrack(vtx[2], eta, thirdDim);
+
+          if (TMath::Abs(vtx[2]) < 10)
+          {
+            fPartEta[0]->Fill(eta);
+
+            if (vtx[2] < 0)
+              fPartEta[1]->Fill(eta);
+            else
+              fPartEta[2]->Fill(eta);
+          }
+          
+          if (TMath::Abs(eta) < 0.5)
+            part05++;
+          if (TMath::Abs(eta) < 1.0) // in the INEL>0, MB WG definition, this is in any case equivalent to <0.8, due to the EsdTrackCuts
+            part10++;
+        }
+        
+        Int_t multAxis = inputCount;
+        if (fMultAxisEta1)
+          multAxis = part10;
+
+        fMultVtx->Fill(multAxis);
+        //if (TMath::Abs(vtx[2]) < 10)
+        //  fMultVtx->Fill(part05);
+
+        // for event count per vertex
+        fdNdEtaAnalysisESD->FillEvent(vtx[2], multAxis);
+
+        // control hist
+        if (multAxis > 0)
+          fEvents->Fill(vtx[2]);
+
+        if (fReadMC)
+        {
+          // from tracks is only done for triggered and vertex reconstructed events
+          for (Int_t i=0; i<inputCount; ++i)
+          {
+            Int_t label = labelArr[i];
+
+            if (label < 0)
+              continue;
+
+            //Printf("Getting particle of track %d", label);
+            TParticle* particle = stack->Particle(label);
+
+            if (!particle)
+            {
+              AliDebug(AliLog::kError, Form("ERROR: Could not retrieve particle %d.", label));
+              continue;
+            }
+
+            Float_t thirdDim = -1;
+            if (fAnalysisMode & AliPWG0Helper::kSPD)
+            {
+              if (fFillPhi)
+              {
+                thirdDim = particle->Phi();
+              }
+              else
+                thirdDim = multAxis;
+            }
+            else
+              thirdDim = particle->Pt();
+
+            Float_t eta = particle->Eta();
+            if (fSymmetrize)
+              eta = TMath::Abs(eta);
+            fdNdEtaAnalysisTracks->FillTrack(vtxMC[2], eta, thirdDim);
+          } // end of track loop
+
+          // for event count per vertex
+          fdNdEtaAnalysisTracks->FillEvent(vtxMC[2], multAxis);
+        }
+      }
     }
+    if (etaArr)
+      delete[] etaArr;
+    if (labelArr)
+      delete[] labelArr;
+    if (thirdDimArr)
+      delete[] thirdDimArr;
   }
 
   if (fReadMC)   // Processing of MC information (optional)
@@ -423,86 +1094,168 @@ void AlidNdEtaTask::Exec(Option_t*)
 
     // get the MC vertex
     AliGenEventHeader* genHeader = header->GenEventHeader();
+    if (!genHeader)
+    {
+      AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
+      return;
+    }
+
     TArrayF vtxMC(3);
     genHeader->PrimaryVertex(vtxMC);
+    
+    // get process type
+    Int_t processType = AliPWG0Helper::GetEventProcessType(fESD, header, stack, fDiffTreatment);
+    AliDebug(AliLog::kDebug+1, Form("Found process type %d", processType));
+
+    if (processType == AliPWG0Helper::kInvalidProcess)
+      AliDebug(AliLog::kError, Form("Unknown process type %d.", processType));
 
     // loop over mc particles
     Int_t nPrim  = stack->GetNprimary();
 
     Int_t nAcceptedParticles = 0;
+    Bool_t oneParticleEvent = kFALSE;
+
+    Int_t nMCparticlesinRange = 0; // number of true particles in my range of interest
 
+    // count particles first, then fill
     for (Int_t iMc = 0; iMc < nPrim; ++iMc)
     {
+      if (!stack->IsPhysicalPrimary(iMc))
+        continue;
+
       //Printf("Getting particle %d", iMc);
       TParticle* particle = stack->Particle(iMc);
 
       if (!particle)
         continue;
 
-      if (AliPWG0Helper::IsPrimaryCharged(particle, nPrim) == kFALSE)
+      if (particle->GetPDG()->Charge() == 0)
         continue;
 
-      AliDebug(AliLog::kDebug+1, Form("Accepted primary %d, unique ID: %d", iMc, particle->GetUniqueID()));
+      //AliDebug(AliLog::kDebug+1, Form("Accepted primary %d, unique ID: %d", iMc, particle->GetUniqueID()));
       Float_t eta = particle->Eta();
-      Float_t pt = particle->Pt();
 
-       // make a rough eta cut (so that nAcceptedParticles is not too far off the true measured value (NB: this histograms are only gathered for comparison))
-      if (TMath::Abs(eta) < 1.5 && pt > 0.3)
-        nAcceptedParticles++;
+      AliDebug(2,Form("particle %d: eta = %f, pt = %f",iMc,particle->Eta(),particle->Pt()));
+
+      // INEL>0: choose one
+      // 1.HL
+      // if (TMath::Abs(eta) < 1.0){
+      //   oneParticleEvent = kTRUE;
+      //   nMCparticlesinRange++;
+      //}
+      // 2.MB Working Group definition
+      if (TMath::Abs(eta) < 0.8 && particle->Pt() > fPtMin){
+       oneParticleEvent = kTRUE;
+       nMCparticlesinRange++;
+      }
 
-      fdNdEtaAnalysis->FillTrack(vtxMC[2], eta, pt);
-      fVertex->Fill(particle->Vx(), particle->Vy(), particle->Vz());
+      // make a rough eta cut (so that nAcceptedParticles is not too far off the true measured value (NB: these histograms are only gathered for comparison))
+      if (TMath::Abs(eta) < 1.5) // && pt > 0.3)
+        nAcceptedParticles++;
+    }
 
-      if (eventTriggered)
-      {
-        fdNdEtaAnalysisTr->FillTrack(vtxMC[2], eta, pt);
-        if (vtxESD)
-          fdNdEtaAnalysisTrVtx->FillTrack(vtxMC[2], eta, pt);
+    if (TMath::Abs(vtxMC[2]) < 10.){
+      // if MC vertex is inside 10
+      if (vtxESD){
+       fVtxEffNum->Fill(nMCparticlesinRange);
       }
-
-      if (TMath::Abs(eta) < 0.8)
-        fPartPt->Fill(pt);
+      fVtxEffDen->Fill(nMCparticlesinRange);
+      if (eventTriggered){
+       if (vtxESD){
+          fVtxTrigEffNum->Fill(nMCparticlesinRange);
+       }
+       fVtxTrigEffDen->Fill(nMCparticlesinRange);
+       fTrigEffNum->Fill(nMCparticlesinRange);
+      }
+      fTrigEffDen->Fill(nMCparticlesinRange);
     }
 
-    fdNdEtaAnalysis->FillEvent(vtxMC[2], nAcceptedParticles);
-    if (eventTriggered)
-    {
-      fdNdEtaAnalysisTr->FillEvent(vtxMC[2], nAcceptedParticles);
-      if (vtxESD)
-        fdNdEtaAnalysisTrVtx->FillEvent(vtxMC[2], nAcceptedParticles);
+    if (oneParticleEvent){
+           fHistEventsMC->Fill(1.1);
     }
+    else{
+           fHistEventsMC->Fill(0.1);
+    }  
 
-    if (eventTriggered && vtxESD)
+    for (Int_t iMc = 0; iMc < nPrim; ++iMc)
     {
-      // from tracks is only done for triggered and vertex reconstructed events
+      if (!stack->IsPhysicalPrimary(iMc))
+        continue;
 
-      for (Int_t i=0; i<inputCount; ++i)
-      {
-        Int_t label = labelArr[i];
+      //Printf("Getting particle %d", iMc);
+      TParticle* particle = stack->Particle(iMc);
 
-        if (label < 0)
-          continue;
+      if (!particle)
+        continue;
+
+      if (particle->GetPDG()->Charge() == 0)
+        continue;
 
-        //Printf("Getting particle of track %d", label);
-        TParticle* particle = stack->Particle(label);
+      Float_t eta = particle->Eta();
+      if (fSymmetrize)
+        eta = TMath::Abs(eta);
+
+      Float_t thirdDim = -1;
 
-        if (!particle)
+      if (fAnalysisMode & AliPWG0Helper::kSPD)
+      {
+        if (fFillPhi)
         {
-          AliDebug(AliLog::kError, Form("ERROR: Could not retrieve particle %d.", label));
-          continue;
+          thirdDim = particle->Phi();
         }
+        else
+          thirdDim = nAcceptedParticles;
+      }
+      else
+        thirdDim = particle->Pt();
+
+      fdNdEtaAnalysis->FillTrack(vtxMC[2], eta, thirdDim);
+
+      if (processType != AliPWG0Helper::kSD)
+        fdNdEtaAnalysisNSD->FillTrack(vtxMC[2], eta, thirdDim);
+
+      if (processType == AliPWG0Helper::kND)
+        fdNdEtaAnalysisND->FillTrack(vtxMC[2], eta, thirdDim);
+      
+      if (oneParticleEvent){
+             AliDebug(3,Form("filling dNdEtaAnalysis object:: vtx = %f, eta = %f, pt = %f",vtxMC[2], eta, thirdDim));
+             fdNdEtaAnalysisOnePart->FillTrack(vtxMC[2], eta, thirdDim);
+      }
+
+      if (eventTriggered)
+      {
+        fdNdEtaAnalysisTr->FillTrack(vtxMC[2], eta, thirdDim);
+        if (vtxESD)
+          fdNdEtaAnalysisTrVtx->FillTrack(vtxMC[2], eta, thirdDim);
+      }
 
-        fdNdEtaAnalysisTracks->FillTrack(vtxMC[2], particle->Eta(), particle->Pt());
-      } // end of track loop
+      if (TMath::Abs(eta) < 1.0 && particle->Pt() > 0 && particle->P() > 0)
+      {
+        //Float_t value = 1. / TMath::TwoPi() / particle->Pt() * particle->Energy() / particle->P();
+        Float_t value = 1;
+        fPartPt->Fill(particle->Pt(), value);
+       if (TMath::Abs(eta) < 0.8 && particle->Pt() > fPtMin){
+         fEtaMC->Fill(eta);
+       }
+      }
+    }
 
-      // for event count per vertex
-      fdNdEtaAnalysisTracks->FillEvent(vtxMC[2], inputCount);
+    fdNdEtaAnalysis->FillEvent(vtxMC[2], nAcceptedParticles);
+    if (processType != AliPWG0Helper::kSD)
+      fdNdEtaAnalysisNSD->FillEvent(vtxMC[2], nAcceptedParticles);
+    if (processType == AliPWG0Helper::kND)
+      fdNdEtaAnalysisND->FillEvent(vtxMC[2], nAcceptedParticles);
+    if (oneParticleEvent)
+      fdNdEtaAnalysisOnePart->FillEvent(vtxMC[2], nAcceptedParticles);
+
+    if (eventTriggered)
+    {
+      fdNdEtaAnalysisTr->FillEvent(vtxMC[2], nAcceptedParticles);
+      if (vtxESD)
+        fdNdEtaAnalysisTrVtx->FillEvent(vtxMC[2], nAcceptedParticles);
     }
   }
-
-  delete[] etaArr;
-  delete[] labelArr;
-  delete[] ptArr;
 }
 
 void AlidNdEtaTask::Terminate(Option_t *)
@@ -511,104 +1264,253 @@ void AlidNdEtaTask::Terminate(Option_t *)
   // a query. It always runs on the client, it can be used to present
   // the results graphically or save the results to file.
 
-  fOutput = dynamic_cast<TList*> (GetOutputData(0));
-  if (!fOutput) {
+  fOutput = dynamic_cast<TList*> (GetOutputData(1));
+  if (!fOutput)
     Printf("ERROR: fOutput not available");
-    return;
-  }
 
-  fdNdEtaAnalysisESD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("fdNdEtaAnalysisESD"));
-  fMult = dynamic_cast<TH1F*> (fOutput->FindObject("fMult"));
-  fMultVtx = dynamic_cast<TH1F*> (fOutput->FindObject("fMultVtx"));
-  fEvents = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_vertex"));
-  fVertexResolution = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_vertex_resolution_z"));
-  for (Int_t i=0; i<3; ++i)
-    fPartEta[i] = dynamic_cast<TH1F*> (fOutput->FindObject(Form("dndeta_check_%d", i)));
+  TH1D* dNdEta = new TH1D("dNdEta","#eta;counts",80, -4, 4);
+  TH1D* dNdEtaMC = new TH1D("dNdEtaMC","#eta,MC;counts",80, -4, 4);
+  if (fOutput)
+  {
+    fdNdEtaAnalysisESD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("fdNdEtaAnalysisESD"));
+    fMult = dynamic_cast<TH1F*> (fOutput->FindObject("fMult"));
+    fMultVtx = dynamic_cast<TH1F*> (fOutput->FindObject("fMultVtx"));
+    for (Int_t i=0; i<3; ++i)
+      fPartEta[i] = dynamic_cast<TH1F*> (fOutput->FindObject(Form("dndeta_check_%d", i)));
+    fEvents = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_vertex"));
+    fVertexResolution = dynamic_cast<TH2F*> (fOutput->FindObject("dndeta_vertex_resolution_z"));
+
+    fVertex = dynamic_cast<TH3F*> (fOutput->FindObject("vertex_check"));
+    fVertexVsMult = dynamic_cast<TH3F*> (fOutput->FindObject("fVertexVsMult"));
+    fPhi = dynamic_cast<TH1F*> (fOutput->FindObject("fPhi"));
+    fRawPt = dynamic_cast<TH1F*> (fOutput->FindObject("fRawPt"));
+    fEtaPhi = dynamic_cast<TH2F*> (fOutput->FindObject("fEtaPhi"));
+    for (Int_t i=0; i<2; ++i)
+      fZPhi[i] = dynamic_cast<TH2F*> (fOutput->FindObject(Form("fZPhi_%d", i)));
+    fModuleMap = dynamic_cast<TH1F*> (fOutput->FindObject("fModuleMap"));
+    fDeltaPhi = dynamic_cast<TH1F*> (fOutput->FindObject("fDeltaPhi"));
+    fDeltaTheta = dynamic_cast<TH1F*> (fOutput->FindObject("fDeltaTheta"));
+    fFiredChips = dynamic_cast<TH2F*> (fOutput->FindObject("fFiredChips"));
+    fTrackletsVsClusters = dynamic_cast<TH2F*> (fOutput->FindObject("fTrackletsVsClusters"));
+    fTrackletsVsUnassigned = dynamic_cast<TH2F*> (fOutput->FindObject("fTrackletsVsUnassigned"));
+    fStats = dynamic_cast<TH1F*> (fOutput->FindObject("fStats"));
+    fStats2 = dynamic_cast<TH2F*> (fOutput->FindObject("fStats2"));
+
+    fEsdTrackCuts = dynamic_cast<AliESDtrackCuts*> (fOutput->FindObject("fEsdTrackCuts"));
+
+    fEta = dynamic_cast<TH1D*>(fOutput->FindObject("fEta"));
+    fEtaMC = dynamic_cast<TH1D*>(fOutput->FindObject("fEtaMC"));
+    fHistEvents = dynamic_cast<TH1D*>(fOutput->FindObject("fHistEvents"));
+    fHistEventsMC = dynamic_cast<TH1D*>(fOutput->FindObject("fHistEventsMC"));
+    fVtxEffDen = dynamic_cast<TH1D*>(fOutput->FindObject("fVtxEffDen"));
+    fVtxEffNum = dynamic_cast<TH1D*>(fOutput->FindObject("fVtxEffNum"));
+    fTrigEffDen = dynamic_cast<TH1D*>(fOutput->FindObject("fTrigEffDen"));
+    fTrigEffNum = dynamic_cast<TH1D*>(fOutput->FindObject("fTrigEffNum"));
+    fVtxTrigEffDen = dynamic_cast<TH1D*>(fOutput->FindObject("fVtxTrigEffDen"));
+    fVtxTrigEffNum = dynamic_cast<TH1D*>(fOutput->FindObject("fVtxTrigEffNum"));
+    Printf("Events selected = %f", fHistEvents->GetBinContent(fHistEvents->GetXaxis()->FindBin(1.1)));
+    Printf("Events selected frm MC = %f", fHistEventsMC->GetBinContent(fHistEventsMC->GetXaxis()->FindBin(1.1)));
+    Float_t nevents = fHistEvents->GetBinContent(fHistEvents->GetXaxis()->FindBin(1.1));
+    Float_t neventsMC = fHistEventsMC->GetBinContent(fHistEventsMC->GetXaxis()->FindBin(1.1));
+    for (Int_t ibin = 1; ibin <= fEta->GetNbinsX(); ibin++){
+      Float_t eta = fEta->GetBinContent(ibin)/nevents/fEta->GetBinWidth(ibin);
+      Float_t etaerr = fEta->GetBinError(ibin)/nevents/fEta->GetBinWidth(ibin);
+      Float_t etaMC = fEtaMC->GetBinContent(ibin)/neventsMC/fEtaMC->GetBinWidth(ibin);
+      Float_t etaerrMC = fEtaMC->GetBinError(ibin)/neventsMC/fEtaMC->GetBinWidth(ibin);
+      dNdEta->SetBinContent(ibin,eta);
+      dNdEta->SetBinError(ibin,etaerr);
+      dNdEtaMC->SetBinContent(ibin,etaMC);
+      dNdEtaMC->SetBinError(ibin,etaerrMC);
+    }
+    new TCanvas("eta", " eta ",50, 50, 550, 550) ;
+    fEta->Draw();
+    new TCanvas("etaMC", " etaMC ",50, 50, 550, 550) ;
+    fEtaMC->Draw();
+    new TCanvas("dNdEta", "#eta;dNdEta ",50, 50, 550, 550) ;
+    dNdEta->Draw();
+    new TCanvas("dNdEtaMC", "#eta,MC;dNdEta ",50, 50, 550, 550) ;
+    dNdEtaMC->Draw();
+    new TCanvas("Events", "Events;Events ",50, 50, 550, 550) ;
+    fHistEvents->Draw();
+    new TCanvas("Events, MC", "Events, MC;Events ",50, 50, 550, 550) ;
+    fHistEventsMC->Draw();
+    
+    TFile* outputFileCheck = new TFile("histogramsCheck.root", "RECREATE");
+    fEta->Write();
+    fEtaMC->Write();
+    dNdEta->Write();
+    dNdEtaMC->Write();
+    outputFileCheck->Write();
+    outputFileCheck->Close();
+
+  }
 
   if (!fdNdEtaAnalysisESD)
   {
     AliDebug(AliLog::kError, "ERROR: fdNdEtaAnalysisESD not available");
-    return;
   }
-
-  if (fMult && fMultVtx)
+  else
   {
-    new TCanvas;
-    fMult->Draw();
-    fMultVtx->SetLineColor(2);
-    fMultVtx->Draw("SAME");
-  }
+    if (fMult && fMultVtx)
+    {
+      new TCanvas;
+      fMult->Draw();
+      fMultVtx->SetLineColor(2);
+      fMultVtx->Draw("SAME");
+    }
 
-  if (fPartEta[0])
-  {
-    Int_t events1 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(-19.9), fEvents->GetXaxis()->FindBin(-0.001));
-    Int_t events2 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(0.001), fEvents->GetXaxis()->FindBin(19.9));
+    if (fFiredChips)
+    {
+      new TCanvas;
+      fFiredChips->Draw("COLZ");
+    }
+
+    if (fPartEta[0])
+    {
+      Int_t events1 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(-9.9), fEvents->GetXaxis()->FindBin(-0.001));
+      Int_t events2 = (Int_t) fEvents->Integral(fEvents->GetXaxis()->FindBin(0.001), fEvents->GetXaxis()->FindBin(9.9));
+
+      Printf("%d events with vertex used", events1 + events2);
+      Printf("%d total events with vertex, %d total triggered events, %d triggered events with 0 multiplicity", (Int_t) fMultVtx->Integral(), (Int_t) fMult->Integral(), (Int_t) fMult->GetBinContent(1));
 
-    Printf("%d events with vertex used", events1 + events2);
+      if (events1 > 0 && events2 > 0)
+      {
+        fPartEta[0]->Scale(1.0 / (events1 + events2));
+        fPartEta[1]->Scale(1.0 / events1);
+        fPartEta[2]->Scale(1.0 / events2);
+
+        for (Int_t i=0; i<3; ++i)
+          fPartEta[i]->Scale(1.0 / fPartEta[i]->GetBinWidth(1));
+
+        new TCanvas("control", "control", 500, 500);
+        for (Int_t i=0; i<3; ++i)
+        {
+          fPartEta[i]->SetLineColor(i+1);
+          fPartEta[i]->Draw((i==0) ? "" : "SAME");
+        }
+      }
+    }
 
-    if (events1 > 0 && events2 > 0)
+    if (fEvents)
+    {
+        new TCanvas("control3", "control3", 500, 500);
+        fEvents->Draw();
+    }
+    
+    if (fStats2)
     {
-    fPartEta[0]->Scale(1.0 / (events1 + events2));
-    fPartEta[1]->Scale(1.0 / events1);
-    fPartEta[2]->Scale(1.0 / events2);
+      new TCanvas;
+      fStats2->Draw("TEXT");
+    }
 
-    for (Int_t i=0; i<3; ++i)
-      fPartEta[i]->Scale(1.0 / fPartEta[i]->GetBinWidth(1));
+    TFile* fout = new TFile("analysis_esd_raw.root", "RECREATE");
+
+    if (fdNdEtaAnalysisESD)
+      fdNdEtaAnalysisESD->SaveHistograms();
+
+    if (fEsdTrackCuts)
+      fEsdTrackCuts->SaveHistograms("esd_track_cuts");
+
+    if (fMult)
+      fMult->Write();
+
+    if (fMultVtx)
+      fMultVtx->Write();
 
-    new TCanvas("control", "control", 500, 500);
     for (Int_t i=0; i<3; ++i)
+      if (fPartEta[i])
+        fPartEta[i]->Write();
+
+    if (fEvents)
+      fEvents->Write();
+
+    if (fVertexResolution)
     {
-      fPartEta[i]->SetLineColor(i+1);
-      fPartEta[i]->Draw((i==0) ? "" : "SAME");
+      fVertexResolution->Write();
+      fVertexResolution->ProjectionX();
+      fVertexResolution->ProjectionY();
     }
-    }
-  }
 
-  if (fEvents)
-  {
-    new TCanvas("control3", "control3", 500, 500);
-    fEvents->Draw();
-  }
+    if (fDeltaPhi)
+      fDeltaPhi->Write();
 
-  TFile* fout = new TFile("analysis_esd_raw.root", "RECREATE");
+    if (fDeltaTheta)
+      fDeltaTheta->Write();
+    
+    if (fPhi)
+      fPhi->Write();
 
-  if (fdNdEtaAnalysisESD)
-    fdNdEtaAnalysisESD->SaveHistograms();
+    if (fRawPt)
+      fRawPt->Write();
 
-  if (fEsdTrackCuts)
-    fEsdTrackCuts->SaveHistograms("esd_tracks_cuts");
+    if (fEtaPhi)
+      fEtaPhi->Write();
 
-  if (fMult)
-    fMult->Write();
+    for (Int_t i=0; i<2; ++i)
+      if (fZPhi[i])
+        fZPhi[i]->Write();
+    
+    if (fModuleMap)
+      fModuleMap->Write();
+    
+    if (fFiredChips)
+      fFiredChips->Write();
 
-  if (fMultVtx)
-    fMultVtx->Write();
+    if (fTrackletsVsClusters)
+      fTrackletsVsClusters->Write();
+    
+    if (fTrackletsVsUnassigned)
+      fTrackletsVsUnassigned->Write();
+    
+    if (fStats)
+      fStats->Write();
 
-  for (Int_t i=0; i<3; ++i)
-    if (fPartEta[i])
-      fPartEta[i]->Write();
+    if (fStats2)
+      fStats2->Write();
+    
+    if (fVertex)
+      fVertex->Write();
 
-  if (fEvents)
-    fEvents->Write();
+    if (fVertexVsMult)
+      fVertexVsMult->Write();
+    
+    if (fVtxEffDen) fVtxEffDen->Write();
+    else Printf("fVtxEffDen not found");
+
+    if (fVtxEffNum) fVtxEffNum->Write();
+    else Printf("fVtxEffNum not found");
+
+    if (fVtxTrigEffNum) fVtxTrigEffNum->Write();
+    else Printf("fVtxTrigEffNum not found");
 
-  if (fVertexResolution)
-    fVertexResolution->Write();
+    if (fVtxTrigEffDen) fVtxTrigEffDen->Write();
+    else Printf("fVtxTrigEffDen not found");
 
-  if (fDeltaPhi)
-    fDeltaPhi->Write();
+    if (fTrigEffNum) fTrigEffNum->Write();
+    else Printf("fTrigEffNum not found");
 
-  fout->Write();
-  fout->Close();
+    if (fTrigEffDen) fTrigEffDen->Write();
+    else Printf("fTrigEffDen not found");
 
-  Printf("Writting result to analysis_esd_raw.root");
+    fout->Write();
+    fout->Close();
+
+    Printf("Writting result to analysis_esd_raw.root");
+  }
 
   if (fReadMC)
   {
-    fdNdEtaAnalysis = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndeta"));
-    fdNdEtaAnalysisTr = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTr"));
-    fdNdEtaAnalysisTrVtx = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTrVtx"));
-    fdNdEtaAnalysisTracks = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTracks"));
-    fPartPt = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_pt"));
+    if (fOutput)
+    {
+      fdNdEtaAnalysis = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndeta"));
+      fdNdEtaAnalysisND = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaND"));
+      fdNdEtaAnalysisNSD = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaNSD"));
+      fdNdEtaAnalysisOnePart = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaOnePart"));
+      fdNdEtaAnalysisTr = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTr"));
+      fdNdEtaAnalysisTrVtx = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTrVtx"));
+      fdNdEtaAnalysisTracks = dynamic_cast<dNdEtaAnalysis*> (fOutput->FindObject("dndetaTracks"));
+      fPartPt = dynamic_cast<TH1F*> (fOutput->FindObject("dndeta_check_pt"));
+    }
 
     if (!fdNdEtaAnalysis || !fdNdEtaAnalysisTr || !fdNdEtaAnalysisTrVtx || !fPartPt)
     {
@@ -617,6 +1519,9 @@ void AlidNdEtaTask::Terminate(Option_t *)
     }
 
     fdNdEtaAnalysis->Finish(0, -1, AlidNdEtaCorrection::kNone);
+    fdNdEtaAnalysisND->Finish(0, -1, AlidNdEtaCorrection::kNone);
+    fdNdEtaAnalysisNSD->Finish(0, -1, AlidNdEtaCorrection::kNone);
+    fdNdEtaAnalysisOnePart->Finish(0, -1, AlidNdEtaCorrection::kNone);
     fdNdEtaAnalysisTr->Finish(0, -1, AlidNdEtaCorrection::kNone);
     fdNdEtaAnalysisTrVtx->Finish(0, -1, AlidNdEtaCorrection::kNone);
     fdNdEtaAnalysisTracks->Finish(0, -1, AlidNdEtaCorrection::kNone);
@@ -628,10 +1533,15 @@ void AlidNdEtaTask::Terminate(Option_t *)
     TFile* fout = new TFile("analysis_mc.root","RECREATE");
 
     fdNdEtaAnalysis->SaveHistograms();
+    fdNdEtaAnalysisND->SaveHistograms();
+    fdNdEtaAnalysisNSD->SaveHistograms();
+    fdNdEtaAnalysisOnePart->SaveHistograms();
     fdNdEtaAnalysisTr->SaveHistograms();
     fdNdEtaAnalysisTrVtx->SaveHistograms();
     fdNdEtaAnalysisTracks->SaveHistograms();
-    fPartPt->Write();
+
+    if (fPartPt)
+      fPartPt->Write();
 
     fout->Write();
     fout->Close();