]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGLF/FORWARD/analysis2/AliCentralMCCorrectionsTask.cxx
Changes for Root6: removing obsolete TH1 functions, corrected EINCLUDE, additional...
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliCentralMCCorrectionsTask.cxx
index 493d99e2e4f9c25241201ddd91fca6ffdb7a67b7..39c193fff3e04d81db108b3333250ab624af82a4 100644 (file)
@@ -12,6 +12,7 @@
 // Corrections used 
 // 
 #include "AliCentralMCCorrectionsTask.h"
+#include "AliCentralCorrectionManager.h"
 #include "AliTriggerAnalysis.h"
 #include "AliPhysicsSelection.h"
 #include "AliLog.h"
@@ -25,6 +26,9 @@
 #include "AliMCEvent.h"
 #include "AliAODForwardMult.h"
 #include "AliCentralCorrSecondaryMap.h"
+#include "AliCentralCorrAcceptance.h"
+#include "AliForwardUtil.h"
+#include "AliMultiplicity.h"
 #include <TH1.h>
 #include <TH2D.h>
 #include <TDirectory.h>
@@ -38,37 +42,18 @@ namespace {
   {
     return Form("nEvents%s%s", (tr ? "Tr" : ""), (vtx ? "Vtx" : ""));
   }
-  /*const char* GetHitsName(UShort_t d, Char_t r) 
-  {
-    return Form("hitsSPD%d%c", d, r);
-  }
-  const char* GetStripsName(UShort_t d, Char_t r)
-  {
-    return Form("stripsSPD%d%c", d, r);
-  }
-  const char* GetPrimaryName(Char_t r, Bool_t trVtx)
-  {
-    return Form("primaries%s%s", 
-               ((r == 'I' || r == 'i') ? "Inner" : "Outer"), 
-               (trVtx ? "TrVtx" : "All"));
-               }
-  */
 }
 
 //====================================================================
 AliCentralMCCorrectionsTask::AliCentralMCCorrectionsTask()
-  : AliAnalysisTaskSE(),
-    fInspector(),
+  : AliBaseMCCorrectionsTask(),
     fTrackDensity(),
-    fVtxBins(0),
-    fFirstEvent(true),
-    fHEvents(0), 
-    fHEventsTr(0), 
-    fHEventsTrVtx(0),
-    fVtxAxis(),
-    fEtaAxis(),
-    fList(),
-    fNPhiBins(20)    
+    fSecCorr(0), 
+    fAccCorr(0),
+    fNPhiBins(20),
+    fEffectiveCorr(true),
+    fEtaCut(1.9),
+    fCorrCut(0.6)
 {
   // 
   // Constructor 
@@ -76,22 +61,19 @@ AliCentralMCCorrectionsTask::AliCentralMCCorrectionsTask()
   // Parameters:
   //    name Name of task 
   //
+  DGUARD(fDebug, 3,"Default CTOR of AliCentralMCCorrectionsTask");
 }
 
 //____________________________________________________________________
 AliCentralMCCorrectionsTask::AliCentralMCCorrectionsTask(const char* name)
-  : AliAnalysisTaskSE(name),
-    fInspector("eventInspector"), 
+  : AliBaseMCCorrectionsTask(name, &(AliCentralCorrectionManager::Instance())),
     fTrackDensity("trackDensity"),
-    fVtxBins(0),
-    fFirstEvent(true),
-    fHEvents(0), 
-    fHEventsTr(0), 
-    fHEventsTrVtx(0),
-    fVtxAxis(10,-10,10), 
-    fEtaAxis(200,-4,6),
-    fList(),
-    fNPhiBins(20)    
+    fSecCorr(0), 
+    fAccCorr(0),
+    fNPhiBins(20),
+    fEffectiveCorr(true),
+    fEtaCut(1.9),
+    fCorrCut(0.6)
 {
   // 
   // Constructor 
@@ -99,397 +81,94 @@ AliCentralMCCorrectionsTask::AliCentralMCCorrectionsTask(const char* name)
   // Parameters:
   //    name Name of task 
   //
-  DefineOutput(1, TList::Class());
-  DefineOutput(2, TList::Class());
-}
-
-//____________________________________________________________________
-AliCentralMCCorrectionsTask::AliCentralMCCorrectionsTask(const AliCentralMCCorrectionsTask& o)
-  : AliAnalysisTaskSE(o),
-    fInspector(o.fInspector),
-    fTrackDensity(),
-    fVtxBins(0),
-    fFirstEvent(o.fFirstEvent),
-    fHEvents(o.fHEvents), 
-    fHEventsTr(o.fHEventsTr), 
-    fHEventsTrVtx(o.fHEventsTrVtx),
-    fVtxAxis(10,-10,10), 
-    fEtaAxis(200,-4,6),
-    fList(o.fList),
-    fNPhiBins(o.fNPhiBins)    
-{
-  // 
-  // Copy constructor 
-  // 
-  // Parameters:
-  //    o Object to copy from 
-  //
-}
-
-//____________________________________________________________________
-AliCentralMCCorrectionsTask&
-AliCentralMCCorrectionsTask::operator=(const AliCentralMCCorrectionsTask& o)
-{
-  // 
-  // Assignment operator 
-  // 
-  // Parameters:
-  //    o Object to assign from 
-  // 
-  // Return:
-  //    Reference to this object 
-  //
-  if (&o == this) return *this; 
-  fInspector         = o.fInspector;
-  fTrackDensity      = o.fTrackDensity;
-  fVtxBins           = o.fVtxBins;
-  fFirstEvent        = o.fFirstEvent;
-  fHEvents           = o.fHEvents;
-  fHEventsTr         = o.fHEventsTr;
-  fHEventsTrVtx      = o.fHEventsTrVtx;
-  SetVertexAxis(o.fVtxAxis);
-  SetEtaAxis(o.fEtaAxis);
-  fNPhiBins          = o.fNPhiBins;
-
-  return *this;
+  DGUARD(fDebug, 3,"Named CTOR of AliCentralMCCorrectionsTask: %s",name);
 }
 
 //____________________________________________________________________
-void
-AliCentralMCCorrectionsTask::Init()
+AliBaseMCCorrectionsTask::VtxBin*
+AliCentralMCCorrectionsTask::CreateVtxBin(Double_t low, Double_t high)
 {
-  // 
-  // Initialize the task 
-  // 
-  //
+  return new AliCentralMCCorrectionsTask::VtxBin(low,high,fEtaAxis,fNPhiBins);
 }
 
 //____________________________________________________________________
-void
-AliCentralMCCorrectionsTask::SetVertexAxis(Int_t nBin, Double_t min, 
-                                          Double_t max)
-{
-  // 
-  // Set the vertex axis to use
-  // 
-  // Parameters:
-  //    nBins Number of bins
-  //    vzMin Least @f$z@f$ coordinate of interation point
-  //    vzMax Largest @f$z@f$ coordinate of interation point
-  //
-  if (max < min) { 
-    Double_t tmp = min;
-    min          = max;
-    max          = tmp;
-  }
-  if (min < -10) 
-    AliWarning(Form("Minimum vertex %f < -10, make sure you want this",min));
-  if (max > +10) 
-    AliWarning(Form("Minimum vertex %f > +10, make sure you want this",max));
-  fVtxAxis.Set(nBin, min, max);
-}
-//____________________________________________________________________
-void
-AliCentralMCCorrectionsTask::SetVertexAxis(const TAxis& axis)
+Bool_t
+AliCentralMCCorrectionsTask::ProcessESD(const AliESDEvent& esd, 
+                                       const AliMCEvent& mc, 
+                                       AliBaseMCCorrectionsTask::VtxBin& bin,
+                                       Double_t          vz)
 {
-  // 
-  // Set the vertex axis to use
-  // 
-  // Parameters:
-  //    axis Axis
-  //
-  SetVertexAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax());
-}
+  AliCentralMCCorrectionsTask::VtxBin& vb = 
+    static_cast<AliCentralMCCorrectionsTask::VtxBin&>(bin);
 
-//____________________________________________________________________
-void
-AliCentralMCCorrectionsTask::SetEtaAxis(Int_t nBin, Double_t min, Double_t max)
-{
-  // 
-  // Set the eta axis to use
-  // 
-  // Parameters:
-  //    nBins Number of bins
-  //    vzMin Least @f$\eta@f$ 
-  //    vzMax Largest @f$\eta@f$ 
-  //
-  if (max < min) { 
-    Double_t tmp = min;
-    min          = max;
-    max          = tmp;
+  // Now process our input data and store in own ESD object 
+  fTrackDensity.Calculate(mc, vz, *vb.fHits, bin.fPrimary);
+  
+  // Get the ESD object
+  const AliMultiplicity* spdmult = esd.GetMultiplicity();
+  // Count number of tracklets per bin 
+  for(Int_t j = 0; j< spdmult->GetNumberOfTracklets();j++) 
+    vb.fClusters->Fill(spdmult->GetEta(j),spdmult->GetPhi(j));
+  //...and then the unused clusters in layer 1 
+  for(Int_t j = 0; j< spdmult->GetNumberOfSingleClusters();j++) {
+    Double_t eta = -TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.));
+    vb.fClusters->Fill(eta, spdmult->GetPhiSingle(j));
   }
-  if (min < -4) 
-    AliWarning(Form("Minimum eta %f < -4, make sure you want this",min));
-  if (max > +6) 
-    AliWarning(Form("Minimum eta %f > +6, make sure you want this",max));
-  fEtaAxis.Set(nBin, min, max);
+  
+  // Count events  
+  bin.fCounts->Fill(0.5);
+  
+  return true;
 }
 //____________________________________________________________________
 void
-AliCentralMCCorrectionsTask::SetEtaAxis(const TAxis& axis)
+AliCentralMCCorrectionsTask::CreateCorrections(TList* results)
 {
-  // 
-  // Set the eta axis to use
-  // 
-  // Parameters:
-  //    axis Axis
-  //
-  SetEtaAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax());
-}
+  fSecCorr = new AliCentralCorrSecondaryMap;
+  fSecCorr->SetVertexAxis(fVtxAxis);
 
-//____________________________________________________________________
-void
-AliCentralMCCorrectionsTask::DefineBins(TList* l)
-{
-  if (!fVtxBins) fVtxBins = new TObjArray(fVtxAxis.GetNbins(), 1);
-  if (fVtxBins->GetEntries() > 0) return;
+  fAccCorr = new AliCentralCorrAcceptance;
+  fAccCorr->SetVertexAxis(fVtxAxis);
 
-  fVtxBins->SetOwner();
-  for (Int_t i = 1; i <= fVtxAxis.GetNbins(); i++) { 
-    Double_t low  = fVtxAxis.GetBinLowEdge(i);
-    Double_t high = fVtxAxis.GetBinUpEdge(i);
-    VtxBin*  bin  = new VtxBin(low, high, fEtaAxis, fNPhiBins);
-    fVtxBins->AddAt(bin, i);
-    bin->DefineOutput(l);
-  }
+  results->Add(fSecCorr);
+  results->Add(fAccCorr);
 }
 
 //____________________________________________________________________
-void
-AliCentralMCCorrectionsTask::UserCreateOutputObjects()
-{
-  // 
-  // Create output objects 
-  // 
-  //
-  fList = new TList;
-  fList->SetOwner();
-  fList->SetName(Form("%sSums", GetName()));
-
-  DefineBins(fList);
-
-  fHEvents = new TH1I(GetEventName(false,false),
-                     "Number of all events", 
-                     fVtxAxis.GetNbins(), 
-                     fVtxAxis.GetXmin(), 
-                     fVtxAxis.GetXmax());
-  fHEvents->SetXTitle("v_{z} [cm]");
-  fHEvents->SetYTitle("# of events");
-  fHEvents->SetFillColor(kBlue+1);
-  fHEvents->SetFillStyle(3001);
-  fHEvents->SetDirectory(0);
-  fList->Add(fHEvents);
-
-  fHEventsTr = new TH1I(GetEventName(true, false), 
-                       "Number of triggered events",
-                       fVtxAxis.GetNbins(), 
-                       fVtxAxis.GetXmin(), 
-                       fVtxAxis.GetXmax());
-  fHEventsTr->SetXTitle("v_{z} [cm]");
-  fHEventsTr->SetYTitle("# of events");
-  fHEventsTr->SetFillColor(kRed+1);
-  fHEventsTr->SetFillStyle(3001);
-  fHEventsTr->SetDirectory(0);
-  fList->Add(fHEventsTr);
-
-  fHEventsTrVtx = new TH1I(GetEventName(true,true),
-                          "Number of events w/trigger and vertex", 
-                          fVtxAxis.GetNbins(), 
-                          fVtxAxis.GetXmin(), 
-                          fVtxAxis.GetXmax());
-  fHEventsTrVtx->SetXTitle("v_{z} [cm]");
-  fHEventsTrVtx->SetYTitle("# of events");
-  fHEventsTrVtx->SetFillColor(kBlue+1);
-  fHEventsTrVtx->SetFillStyle(3001);
-  fHEventsTrVtx->SetDirectory(0);
-  fList->Add(fHEventsTrVtx);
-
-  // Copy axis objects to output 
-  TH1* vtxAxis = new TH1D("vtxAxis", "Vertex axis", 
-                         fVtxAxis.GetNbins(), 
-                         fVtxAxis.GetXmin(), 
-                         fVtxAxis.GetXmax());
-  TH1* etaAxis = new TH1D("etaAxis", "Eta axis", 
-                         fEtaAxis.GetNbins(), 
-                         fEtaAxis.GetXmin(), 
-                         fEtaAxis.GetXmax());
-  fList->Add(vtxAxis);
-  fList->Add(etaAxis);
-
-  AliInfo(Form("Initialising sub-routines: %p, %p", 
-              &fInspector, &fTrackDensity));
-  fInspector.DefineOutput(fList);
-  fInspector.Init(fVtxAxis);
-  fTrackDensity.DefineOutput(fList);
-
-  PostData(1, fList);
-}
-//____________________________________________________________________
-void
-AliCentralMCCorrectionsTask::UserExec(Option_t*)
+Bool_t 
+AliCentralMCCorrectionsTask::FinalizeVtxBin(AliBaseMCCorrectionsTask::VtxBin* 
+                                           bin,  UShort_t iVz) 
 {
-  // 
-  // Process each event 
-  // 
-  // Parameters:
-  //    option Not used
-  //  
-
-  // Get the input data - MC event
-  AliMCEvent*  mcEvent = MCEvent();
-  if (!mcEvent) { 
-    AliWarning("No MC event found");
-    return;
-  }
-
-  // Get the input data - ESD event
-  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
-  if (!esd) { 
-    AliWarning("No ESD event found for input event");
-    return;
-  }
-
-  //--- Read run information -----------------------------------------
-  if (fFirstEvent && esd->GetESDRun()) {
-    fInspector.ReadRunDetails(esd);
-    
-    AliInfo(Form("Initializing with parameters from the ESD:\n"
-                "         AliESDEvent::GetBeamEnergy()   ->%f\n"
-                "         AliESDEvent::GetBeamType()     ->%s\n"
-                "         AliESDEvent::GetCurrentL3()    ->%f\n"
-                "         AliESDEvent::GetMagneticField()->%f\n"
-                "         AliESDEvent::GetRunNumber()    ->%d\n",
-                esd->GetBeamEnergy(),
-                esd->GetBeamType(),
-                esd->GetCurrentL3(),
-                esd->GetMagneticField(),
-                esd->GetRunNumber()));
-
-    Print();
-    fFirstEvent = false;
-  }
-
-  // Some variables 
-  UInt_t   triggers; // Trigger bits
-  Bool_t   lowFlux;  // Low flux flag
-  UShort_t iVz;      // Vertex bin from ESD
-  Double_t vZ;       // Z coordinate from ESD
-  Double_t cent;     // Centrality 
-  UShort_t iVzMc;    // Vertex bin from MC
-  Double_t vZMc;     // Z coordinate of IP vertex from MC
-  Double_t b;        // Impact parameter
-  Int_t    nPart;    // Number of participants 
-  Int_t    nBin;     // Number of binary collisions 
-  Double_t phiR;     // Reaction plane from MC
-  UShort_t nClusters;// Number of SPD clusters 
-  // Process the data 
-  UInt_t retESD = fInspector.Process(esd, triggers, lowFlux, iVz, vZ, 
-                                    cent, nClusters);
-  fInspector.ProcessMC(mcEvent, triggers, iVzMc, vZMc, 
-                      b, nPart, nBin, phiR);
-
-  Bool_t isInel   = triggers & AliAODForwardMult::kInel;
-  Bool_t hasVtx   = retESD == AliFMDMCEventInspector::kOk;
-
-  // Fill the event count histograms 
-  if (isInel)           fHEventsTr->Fill(vZMc);
-  if (isInel && hasVtx) fHEventsTrVtx->Fill(vZMc);
-  fHEvents->Fill(vZMc);
-
-  // Now find our vertex bin object 
-  VtxBin* bin = 0;
-  if (iVzMc > 0 && iVzMc <= fVtxAxis.GetNbins()) 
-    bin = static_cast<VtxBin*>(fVtxBins->At(iVzMc));
-  if (!bin) { 
-    // AliError(Form("No vertex bin object @ %d (%f)", iVzMc, vZMc));
-    return;
-  }
-
-  // Now process our input data and store in own ESD object 
-  fTrackDensity.Calculate(*mcEvent, vZMc, *bin->fHits, bin->fPrimary);
-  bin->fCounts->Fill(0.5);
   
-  PostData(1, fList);
-}
-
-//____________________________________________________________________
-void
-AliCentralMCCorrectionsTask::Terminate(Option_t*)
-{
-  // 
-  // End of job
-  // 
-  // Parameters:
-  //    option Not used 
-  //
-
-  fList = dynamic_cast<TList*>(GetOutputData(1));
-  if (!fList) {
-    AliError("No output list defined");
-    return;
-  }
-
-  DefineBins(fList);
-
-  // Output list 
-  TList* output = new TList;
-  output->SetOwner();
-  output->SetName(Form("%sResults", GetName()));
-
-  // --- Fill correction object --------------------------------------
-  AliCentralCorrSecondaryMap* corr = new AliCentralCorrSecondaryMap;
-  corr->SetVertexAxis(fVtxAxis);
-  // corr->SetEtaAxis(fEtaAxis);
-
-  TIter     next(fVtxBins);
-  VtxBin*   bin = 0;
-  UShort_t  iVz = 1;
-  while ((bin = static_cast<VtxBin*>(next()))) 
-    bin->Finish(fList, output, iVz++, corr);
-
-  output->Add(corr);
-
-  PostData(2, output);
+  AliCentralMCCorrectionsTask::VtxBin* vb = 
+    static_cast<AliCentralMCCorrectionsTask::VtxBin*>(bin);
+  vb->Terminate(fList, fResults, iVz, fEffectiveCorr, 
+               fEtaCut, fCorrCut, fSecCorr,fAccCorr);
+  return true;
 }
+                                           
 
 //____________________________________________________________________
 void
 AliCentralMCCorrectionsTask::Print(Option_t* option) const
 {
-  std::cout << ClassName() << "\n"
-           << "  Vertex bins:      " << fVtxAxis.GetNbins() << '\n'
-           << "  Vertex range:     [" << fVtxAxis.GetXmin() 
-           << "," << fVtxAxis.GetXmax() << "]\n"
-           << "  Eta bins:         " << fEtaAxis.GetNbins() << '\n'
-           << "  Eta range:        [" << fEtaAxis.GetXmin() 
-           << "," << fEtaAxis.GetXmax() << "]\n"
-           << "  # of phi bins:    " << fNPhiBins 
+  AliBaseMCCorrectionsTask::Print(option);
+  std::cout << "  # of phi bins:    " << fNPhiBins << "\n"
+           << "  Effective corr.:  " << fEffectiveCorr << "\n"
+           << "  Eta cut-off:      " << fEtaCut << "\n"
+           << "  Acceptance cut:   " << fCorrCut 
            << std::endl;
   gROOT->IncreaseDirLevel();
-  fInspector.Print(option);
   fTrackDensity.Print(option);
   gROOT->DecreaseDirLevel();
 }
 
 //====================================================================
-const char*
-AliCentralMCCorrectionsTask::VtxBin::BinName(Double_t low, 
-                                            Double_t high) 
-{
-  static TString buf;
-  buf = Form("vtx%+05.1f_%+05.1f", low, high);
-  buf.ReplaceAll("+", "p");
-  buf.ReplaceAll("-", "m");
-  buf.ReplaceAll(".", "d");
-  return buf.Data();
-}
-
-
-//____________________________________________________________________
 AliCentralMCCorrectionsTask::VtxBin::VtxBin()
-  : fHits(0), 
-    fPrimary(0),
-    fCounts(0)
+  : AliBaseMCCorrectionsTask::VtxBin(),
+    fHits(0), 
+    fClusters(0)
 {
 }
 //____________________________________________________________________
@@ -497,121 +176,164 @@ AliCentralMCCorrectionsTask::VtxBin::VtxBin(Double_t     low,
                                            Double_t     high, 
                                            const TAxis& axis,
                                            UShort_t     nPhi)
-  : TNamed(BinName(low, high), 
-          Form("%+5.1fcm<v_{z}<%+5.1fcm", low, high)),
+  : AliBaseMCCorrectionsTask::VtxBin(low, high, axis, nPhi),
     fHits(0), 
-    fPrimary(0),
-    fCounts(0)
+    fClusters(0)
 {
-  fPrimary = new TH2D("primary", "Primaries", 
-                     axis.GetNbins(), axis.GetXmin(), axis.GetXmax(), 
-                     nPhi, 0, 2*TMath::Pi());
-  fPrimary->SetXTitle("#eta");
-  fPrimary->SetYTitle("#varphi [radians]");
-  fPrimary->Sumw2();
-  fPrimary->SetDirectory(0);
-
   fHits = static_cast<TH2D*>(fPrimary->Clone("hits"));
   fHits->SetTitle("Hits");
   fHits->SetDirectory(0);
 
-  fCounts = new TH1D("counts", "Counts", 1, 0, 1);
-  fCounts->SetXTitle("Events");
-  fCounts->SetYTitle("# of Events");
-  fCounts->SetDirectory(0);
-}
-
-//____________________________________________________________________
-AliCentralMCCorrectionsTask::VtxBin::VtxBin(const VtxBin& o)
-  : TNamed(o),
-    fHits(0),
-    fPrimary(0), 
-    fCounts(0)
-{
-  if (o.fHits) {
-    fHits = static_cast<TH2D*>(o.fHits->Clone());
-    fHits->SetDirectory(0);
-  }
-  if (o.fPrimary) {
-    fPrimary = static_cast<TH2D*>(o.fPrimary->Clone());
-    fPrimary->SetDirectory(0);
-  }
-  if (o.fCounts) {
-    fCounts = static_cast<TH1D*>(o.fCounts->Clone());
-    fCounts->SetDirectory(0);
-  }
+  fClusters = static_cast<TH2D*>(fPrimary->Clone("clusters"));
+  fClusters->SetTitle("Clusters");
+  fClusters->SetDirectory(0);
 }
 
-//____________________________________________________________________
-AliCentralMCCorrectionsTask::VtxBin&
-AliCentralMCCorrectionsTask::VtxBin::operator=(const VtxBin& o)
-{
-  if (&o == this) return *this; 
-  TNamed::operator=(o);
-  fHits    = 0;
-  fPrimary = 0;
-  fCounts  = 0;
-  if (o.fHits) {
-    fHits = static_cast<TH2D*>(o.fHits->Clone());
-    fHits->SetDirectory(0);
-  }
-  if (o.fPrimary) {
-    fPrimary = static_cast<TH2D*>(o.fPrimary->Clone());
-    fPrimary->SetDirectory(0);
-  }
-  if (o.fCounts) {
-    fCounts = static_cast<TH1D*>(o.fCounts->Clone());
-    fCounts->SetDirectory(0);
-  }
-  return *this;
-}
 
 //____________________________________________________________________
-void
-AliCentralMCCorrectionsTask::VtxBin::DefineOutput(TList* l)
+TList*
+AliCentralMCCorrectionsTask::VtxBin::CreateOutputObjects(TList* l)
 {
-  TList* d = new TList;
-  d->SetName(GetName());
-  d->SetOwner();
-  l->Add(d);
+  TList* d = AliBaseMCCorrectionsTask::VtxBin::CreateOutputObjects(l);
 
   d->Add(fHits);
-  d->Add(fPrimary);
-  d->Add(fCounts);
+  d->Add(fClusters);
+
+  return d;
 }
 //____________________________________________________________________
 void
-AliCentralMCCorrectionsTask::VtxBin::Finish(const TList* input, 
-                                           TList* output, 
-                                           UShort_t iVz, 
-                                           AliCentralCorrSecondaryMap* map)
+AliCentralMCCorrectionsTask::VtxBin::Terminate(const TList* input, 
+                                              TList* output, 
+                                              UShort_t iVz, 
+                                              Bool_t effectiveCorr,
+                                              Double_t etaCut,
+                                              Double_t accCut,
+                                              AliCentralCorrSecondaryMap* map,
+                                              AliCentralCorrAcceptance* acorr)
 {
   TList* out = new TList;
   out->SetName(GetName());
   out->SetOwner();
   output->Add(out);
-  
+
   TList* l = static_cast<TList*>(input->FindObject(GetName()));
   if (!l) { 
     AliError(Form("List %s not found in %s", GetName(), input->GetName()));
     return;
   }
 
-
+  // Get the sums 
   TH2D*   hits  = static_cast<TH2D*>(l->FindObject("hits"));
+  TH2D*   clus  = static_cast<TH2D*>(l->FindObject("clusters"));
   TH2D*   prim  = static_cast<TH2D*>(l->FindObject("primary"));
   if (!hits || !prim) {
     AliError(Form("Missing histograms: %p, %p", hits, prim));
     return;
   }
 
-  TH2D* h = static_cast<TH2D*>(hits->Clone("bgCorr"));
-  h->SetDirectory(0);
-  h->Divide(prim);
-
-  map->SetCorrection(iVz, h);
+  // Clone cluster and hit map
+  TH2D* secMapEff = static_cast<TH2D*>(clus->Clone("secMapEff"));
+  TH2D* secMapHit = static_cast<TH2D*>(hits->Clone("secMapHit"));
+  secMapEff->SetTitle("2^{nd} map from clusters");
+  secMapEff->SetDirectory(0);
+  secMapHit->SetTitle("2^{nd} map from MC hits");
+  secMapHit->SetDirectory(0);
+
+  // Divide cluster and hit map with number of primaries 
+  secMapEff->Divide(prim);
+  secMapHit->Divide(prim);
+
+  // Create acceptance histograms 
+  TH1D* accEff = new TH1D("accEff",
+                         "Acceptance correction for SPD (from 2^{nd} map)" ,
+                         fPrimary->GetXaxis()->GetNbins(), 
+                         fPrimary->GetXaxis()->GetXmin(), 
+                         fPrimary->GetXaxis()->GetXmax());
+  TH1D* accHit = static_cast<TH1D*>(accEff->Clone("accHit"));
+  accHit->SetTitle("Acceptance correction for SPD (from clusters)");
+
+  // Diagnostics histogra, 
+  TH2*  dia    = static_cast<TH2D*>(clus->Clone("diagnostics"));
+  dia->SetTitle("Scaled cluster density");
+
+  // Total number of channels along phi and # of eta bins
+  Int_t nTotal = secMapHit->GetNbinsY();
+  Int_t nEta   = secMapHit->GetNbinsX();
+
+  for(Int_t xx = 1; xx <= nEta; xx++) {
+    Double_t eta = secMapHit->GetXaxis()->GetBinCenter(xx);
+    Bool_t   ins = TMath::Abs(eta) <= etaCut;
+    Double_t mm  = 0;
+    if (ins) {
+      // Find the maximum cluster signal in this phi range 
+      for (Int_t yy = 1; yy <= nTotal; yy++) { 
+       Double_t c = clus->GetBinContent(xx,yy);
+       mm         = TMath::Max(mm, c);
+      }
+    }
+    // Count number of phi bins with enough clusters or high enough
+    // correction.
+    Int_t nOKEff    = 0;
+    Int_t nOKHit    = 0;
+    for(Int_t yy = 1; yy <=nTotal; yy++) {
+      if (!ins) { // Not inside Eta cut
+       secMapEff->SetBinContent(xx,yy,0.); 
+       secMapEff->SetBinError(xx,yy,0.); 
+       secMapHit->SetBinContent(xx,yy,0.); 
+       secMapHit->SetBinError(xx,yy,0.); 
+       dia->SetBinContent(xx,yy,0);
+       continue;
+      }
+
+      // Check if the background correction is large enough, or zero map
+      if(secMapEff->GetBinContent(xx,yy) > 0.9) {
+       // acc->Fill(h->GetXaxis()->GetBinCenter(xx));
+       nOKEff++;
+      }
+      else {
+       secMapEff->SetBinContent(xx,yy,0.); 
+       secMapEff->SetBinError(xx,yy,0.); 
+      }
+
+      // Check if the number of cluster is large enough, or zero map
+      Double_t c = clus->GetBinContent(xx,yy);
+      Double_t s = (mm < 1e-6) ? 0 : c / mm;
+      dia->SetBinContent(xx,yy,s);
+      if (s >= accCut) {
+       nOKHit++;
+      }
+      else {
+       secMapHit->SetBinContent(xx,yy,0);
+       secMapHit->SetBinError(xx,yy,0);
+      }
+    }
+
+    // Calculate acceptance as ratio of bins with enough clusters and
+    // total number of phi bins.
+    Double_t accXX = float(nOKHit) / nTotal;
+    if (accXX < 0.2) accXX = 0;
+    accHit->SetBinContent(xx, accXX);
+
+    // Calculate acceptance as ratio of bins with large enough
+    // correction and total number of phi bins.
+    accXX = float(nOKEff) / nTotal;
+    if (accXX < 0.2) accXX = 0;
+    accEff->SetBinContent(xx, accXX);
+  }
 
-  out->Add(h);
+  TH2D* secMap    = (effectiveCorr ? secMapEff : secMapHit);
+  TH2D* secMapAlt = (effectiveCorr ? secMapHit : secMapEff);
+  TH1D* acc       = (effectiveCorr ? accEff    : accHit);
+  TH1D* accAlt    = (effectiveCorr ? accHit    : accEff);
+  out->Add(secMap->Clone("secMap"));
+  out->Add(secMapAlt->Clone());
+  out->Add(acc->Clone("acc"));
+  out->Add(accAlt->Clone());
+  out->Add(dia->Clone());
+
+  map->SetCorrection(iVz, secMap);
+  acorr->SetCorrection(iVz, acc);
 }
 
 //