[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliCentralMCCorrectionsTask.cxx

// 
// Calculate the corrections in the central regions
// 
// Inputs: 
//   - AliESDEvent 
//
// Outputs: 
//   - AliAODCentralMult 
// 
// Histograms 
//   
// Corrections used 
// 
#include "AliCentralMCCorrectionsTask.h"
#include "AliTriggerAnalysis.h"
#include "AliPhysicsSelection.h"
#include "AliLog.h"
#include "AliHeader.h"
#include "AliGenEventHeader.h"
#include "AliESDEvent.h"
#include "AliAODHandler.h"
#include "AliMultiplicity.h"
#include "AliInputEventHandler.h"
#include "AliStack.h"
#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>
#include <TList.h>
#include <TROOT.h>
#include <iostream>

//====================================================================
namespace {
  const char* GetEventName(Bool_t tr, Bool_t vtx) 
  {
    return Form("nEvents%s%s", (tr ? "Tr" : ""), (vtx ? "Vtx" : ""));
  }
}

//====================================================================
AliCentralMCCorrectionsTask::AliCentralMCCorrectionsTask()
  : AliAnalysisTaskSE(),
    fInspector(),
    fTrackDensity(),
    fVtxBins(0),
    fFirstEvent(true),
    fHEvents(0), 
    fHEventsTr(0), 
    fHEventsTrVtx(0),
    fVtxAxis(),
    fEtaAxis(),
    fList(),
    fNPhiBins(20),
    fEffectiveCorr(true)
{
  // 
  // Constructor 
  // 
  // Parameters:
  //    name Name of task 
  //
  DGUARD(fDebug, 3,"Default CTOR of AliCentralMCCorrectionsTask");
}

//____________________________________________________________________
AliCentralMCCorrectionsTask::AliCentralMCCorrectionsTask(const char* name)
  : AliAnalysisTaskSE(name),
    fInspector("eventInspector"), 
    fTrackDensity("trackDensity"),
    fVtxBins(0),
    fFirstEvent(true),
    fHEvents(0), 
    fHEventsTr(0), 
    fHEventsTrVtx(0),
    fVtxAxis(10,-10,10), 
    fEtaAxis(200,-4,6),
    fList(),
    fNPhiBins(20),
    fEffectiveCorr(true)
{
  // 
  // Constructor 
  // 
  // Parameters:
  //    name Name of task 
  //
  DGUARD(fDebug, 3,"Named CTOR of AliCentralMCCorrectionsTask: %s",name);
  DefineOutput(1, TList::Class());
  DefineOutput(2, TList::Class());
  fBranchNames = 
    "ESD:AliESDRun.,AliESDHeader.,AliMultiplicity.,"
    "SPDVertex.,PrimaryVertex.";
}

//____________________________________________________________________
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),
    fEffectiveCorr(o.fEffectiveCorr)
{
  // 
  // Copy constructor 
  // 
  // Parameters:
  //    o Object to copy from 
  //
  DGUARD(fDebug, 3,"Copy CTOR of AliCentralMCCorrectionsTask");
}

//____________________________________________________________________
AliCentralMCCorrectionsTask&
AliCentralMCCorrectionsTask::operator=(const AliCentralMCCorrectionsTask& o)
{
  // 
  // Assignment operator 
  // 
  // Parameters:
  //    o Object to assign from 
  // 
  // Return:
  //    Reference to this object 
  //
  DGUARD(fDebug,3,"Assignment of AliCentralMCCorrectionsTask");
  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;
  fEffectiveCorr     = o.fEffectiveCorr;

  return *this;
}

//____________________________________________________________________
void
AliCentralMCCorrectionsTask::Init()
{
  // 
  // Initialize the task 
  // 
  //
  DGUARD(fDebug,1,"Initialize AliCentralMCCorrectionsTask");
}

//____________________________________________________________________
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
  //
  DGUARD(fDebug,3,"Set vertex axis AliCentralMCCorrectionsTask [%d,%f,%f]",
	 nBin, min, max);
  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)
{
  // 
  // Set the vertex axis to use
  // 
  // Parameters:
  //    axis Axis
  //
  SetVertexAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax());
}

//____________________________________________________________________
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$ 
  //
  DGUARD(fDebug,3,"Set eta axis AliCentralMCCorrectionsTask [%d,%f,%f]",
	 nBin, min, max);
  if (max < min) { 
    Double_t tmp = min;
    min          = max;
    max          = tmp;
  }
  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);
}
//____________________________________________________________________
void
AliCentralMCCorrectionsTask::SetEtaAxis(const TAxis& axis)
{
  // 
  // Set the eta axis to use
  // 
  // Parameters:
  //    axis Axis
  //
  SetEtaAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax());
}

//____________________________________________________________________
void
AliCentralMCCorrectionsTask::DefineBins(TList* l)
{
  DGUARD(fDebug,1,"Define bins in AliCentralMCCorrectionsTask");
  if (!fVtxBins) fVtxBins = new TObjArray(fVtxAxis.GetNbins(), 1);
  if (fVtxBins->GetEntries() > 0) return;

  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->CreateOutputObjects(l);
  }
}

//____________________________________________________________________
void
AliCentralMCCorrectionsTask::UserCreateOutputObjects()
{
  // 
  // Create output objects 
  // 
  //
  DGUARD(fDebug,1,"Create user output for AliCentralMCCorrectionsTask");
  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.CreateOutputObjects(fList);
  fTrackDensity.CreateOutputObjects(fList);

  PostData(1, fList);
}
//____________________________________________________________________
void
AliCentralMCCorrectionsTask::UserExec(Option_t*)
{
  // 
  // Process each event 
  // 
  // Parameters:
  //    option Not used
  //  

  DGUARD(fDebug,1,"AliCentralMCCorrectionsTask process an event");
  // 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()));

    fInspector.SetupForData(fVtxAxis);

    Print();
    fFirstEvent = false;
  }

  // Some variables 
  UInt_t   triggers  = 0;    // Trigger bits
  Bool_t   lowFlux   = true; // Low flux flag
  UShort_t iVz       = 0;    // Vertex bin from ESD
  TVector3 ip;               // Z coordinate from ESD
  Double_t cent      = -1;   // Centrality 
  UShort_t iVzMc     = 0;    // Vertex bin from MC
  Double_t vZMc      = 1000; // Z coordinate of IP vertex from MC
  Double_t b         = -1;   // Impact parameter
  Double_t cMC       = -1;   // Centrality estimate from b
  Int_t    nPart     = -1;   // Number of participants 
  Int_t    nBin      = -1;   // Number of binary collisions 
  Double_t phiR      = 1000; // Reaction plane from MC
  UShort_t nClusters = 0;    // Number of SPD clusters 
  // Process the data 
  UInt_t retESD = fInspector.Process(esd, triggers, lowFlux, iVz, ip, 
				     cent, nClusters);
  fInspector.ProcessMC(mcEvent, triggers, iVzMc, vZMc, 
		       b, cMC, 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);

  // Get the ESD object
  const AliMultiplicity* spdmult = esd->GetMultiplicity();

  // Count number of tracklets per bin 
  for(Int_t j = 0; j< spdmult->GetNumberOfTracklets();j++) 
    bin->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.));
     bin->fClusters->Fill(eta, spdmult->GetPhiSingle(j));
  }

  // Count events  
  bin->fCounts->Fill(0.5);
  
  PostData(1, fList);
}

//____________________________________________________________________
void
AliCentralMCCorrectionsTask::Terminate(Option_t*)
{
  // 
  // End of job
  // 
  // Parameters:
  //    option Not used 
  //
  DGUARD(fDebug,1,"AliCentralMCCorrectionsTask analyse merged output");

  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);

  AliCentralCorrAcceptance* acorr = new AliCentralCorrAcceptance;
  acorr->SetVertexAxis(fVtxAxis);

  TIter     next(fVtxBins);
  VtxBin*   bin = 0;
  UShort_t  iVz = 1;
  while ((bin = static_cast<VtxBin*>(next()))) 
    bin->Terminate(fList, output, iVz++, fEffectiveCorr, corr,acorr);

  output->Add(corr);
  output->Add(acorr);

  PostData(2, output);
}

//____________________________________________________________________
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 
	    << 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), 
    fClusters(0),
    fPrimary(0),
    fCounts(0)
{
}
//____________________________________________________________________
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)),
    fHits(0), 
    fClusters(0),
    fPrimary(0),
    fCounts(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);

  fClusters = static_cast<TH2D*>(fPrimary->Clone("clusters"));
  fClusters->SetTitle("Clusters");
  fClusters->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),
    fClusters(0),
    fPrimary(0), 
    fCounts(0)
{
  if (o.fHits) {
    fHits = static_cast<TH2D*>(o.fHits->Clone());
    fHits->SetDirectory(0);
  }
  if (o.fClusters) {
    fClusters = static_cast<TH2D*>(o.fClusters->Clone());
    fClusters->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);
  }
}

//____________________________________________________________________
AliCentralMCCorrectionsTask::VtxBin&
AliCentralMCCorrectionsTask::VtxBin::operator=(const VtxBin& o)
{
  if (&o == this) return *this; 
  TNamed::operator=(o);
  fHits     = 0;
  fPrimary  = 0;
  fClusters = 0;
  fCounts   = 0;
  if (o.fHits) {
    fHits = static_cast<TH2D*>(o.fHits->Clone());
    fHits->SetDirectory(0);
  }
  if (o.fClusters) {
    fClusters = static_cast<TH2D*>(o.fClusters->Clone());
    fClusters->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::CreateOutputObjects(TList* l)
{
  TList* d = new TList;
  d->SetName(GetName());
  d->SetOwner();
  l->Add(d);

  d->Add(fHits);
  d->Add(fClusters);
  d->Add(fPrimary);
  d->Add(fCounts);
}
//____________________________________________________________________
void
AliCentralMCCorrectionsTask::VtxBin::Terminate(const TList* input, 
					    TList* output, 
					    UShort_t iVz, 
					    Bool_t effectiveCorr,
					    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;
  }


  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 = 0;
  if (effectiveCorr) h = static_cast<TH2D*>(clus->Clone("bgCorr"));
  else               h = static_cast<TH2D*>(hits->Clone("bgCorr"));
  h->SetDirectory(0);
  h->Divide(prim);

  TH1D* acc = new TH1D(Form("SPDacc_vrtbin_%d",iVz),
		       "Acceptance correction for SPD" ,
		       fPrimary->GetXaxis()->GetNbins(), 
		       fPrimary->GetXaxis()->GetXmin(), 
		       fPrimary->GetXaxis()->GetXmax());
  TH1F* accden = static_cast<TH1F*>(acc->Clone(Form("%s_den",acc->GetName())));

  for(Int_t xx = 1; xx <=h->GetNbinsX(); xx++) {
    for(Int_t yy = 1; yy <=h->GetNbinsY(); yy++) {
      if(TMath::Abs(h->GetXaxis()->GetBinCenter(xx)) > 1.9) {
	h->SetBinContent(xx,yy,0.); 
	h->SetBinError(xx,yy,0.); 
      }
      if(h->GetBinContent(xx,yy) > 0.9) 
	acc->Fill(h->GetXaxis()->GetBinCenter(xx));
      else {
	h->SetBinContent(xx,yy,0.); 
	h->SetBinError(xx,yy,0.); 
      }
      accden->Fill(h->GetXaxis()->GetBinCenter(xx));
    }
  }
  acc->Divide(accden);

  map->SetCorrection(iVz, h);
  acorr->SetCorrection(iVz, acc);
  out->Add(h);
  out->Add(acc);
}

//
// EOF
//
Commit	Line	Data
600b313f	1	//
	2	// Calculate the corrections in the central regions
	3	//
	4	// Inputs:
	5	// - AliESDEvent
	6	//
	7	// Outputs:
	8	// - AliAODCentralMult
	9	//
	10	// Histograms
	11	//
	12	// Corrections used
	13	//
	14	#include "AliCentralMCCorrectionsTask.h"
	15	#include "AliTriggerAnalysis.h"
	16	#include "AliPhysicsSelection.h"
	17	#include "AliLog.h"
	18	#include "AliHeader.h"
	19	#include "AliGenEventHeader.h"
	20	#include "AliESDEvent.h"
	21	#include "AliAODHandler.h"
	22	#include "AliMultiplicity.h"
	23	#include "AliInputEventHandler.h"
	24	#include "AliStack.h"
	25	#include "AliMCEvent.h"
	26	#include "AliAODForwardMult.h"
	27	#include "AliCentralCorrSecondaryMap.h"
ea6f96fb	28	#include "AliCentralCorrAcceptance.h"
f53fb4f6	29	#include "AliForwardUtil.h"
ea6f96fb	30	#include "AliMultiplicity.h"
600b313f	31	#include <TH1.h>
	32	#include <TH2D.h>
	33	#include <TDirectory.h>
	34	#include <TList.h>
	35	#include <TROOT.h>
	36	#include <iostream>
	37
	38	//====================================================================
	39	namespace {
	40	const char* GetEventName(Bool_t tr, Bool_t vtx)
	41	{
	42	return Form("nEvents%s%s", (tr ? "Tr" : ""), (vtx ? "Vtx" : ""));
	43	}
600b313f	44	}
	45
	46	//====================================================================
	47	AliCentralMCCorrectionsTask::AliCentralMCCorrectionsTask()
	48	: AliAnalysisTaskSE(),
	49	fInspector(),
	50	fTrackDensity(),
	51	fVtxBins(0),
	52	fFirstEvent(true),
	53	fHEvents(0),
	54	fHEventsTr(0),
	55	fHEventsTrVtx(0),
	56	fVtxAxis(),
	57	fEtaAxis(),
	58	fList(),
ea6f96fb	59	fNPhiBins(20),
ea6f96fb	60	fEffectiveCorr(true)
600b313f	61	{
	62	//
	63	// Constructor
	64	//
	65	// Parameters:
	66	// name Name of task
	67	//
5ca83fee	68	DGUARD(fDebug, 3,"Default CTOR of AliCentralMCCorrectionsTask");
600b313f	69	}
	70
	71	//____________________________________________________________________
	72	AliCentralMCCorrectionsTask::AliCentralMCCorrectionsTask(const char* name)
	73	: AliAnalysisTaskSE(name),
	74	fInspector("eventInspector"),
	75	fTrackDensity("trackDensity"),
	76	fVtxBins(0),
	77	fFirstEvent(true),
	78	fHEvents(0),
	79	fHEventsTr(0),
	80	fHEventsTrVtx(0),
	81	fVtxAxis(10,-10,10),
	82	fEtaAxis(200,-4,6),
	83	fList(),
ea6f96fb	84	fNPhiBins(20),
ea6f96fb	85	fEffectiveCorr(true)
600b313f	86	{
	87	//
	88	// Constructor
	89	//
	90	// Parameters:
	91	// name Name of task
	92	//
5ca83fee	93	DGUARD(fDebug, 3,"Named CTOR of AliCentralMCCorrectionsTask: %s",name);
600b313f	94	DefineOutput(1, TList::Class());
600b313f	95	DefineOutput(2, TList::Class());
3ecc3599	96	fBranchNames =
	97	"ESD:AliESDRun.,AliESDHeader.,AliMultiplicity.,"
	98	"SPDVertex.,PrimaryVertex.";
600b313f	99	}
	100
	101	//____________________________________________________________________
	102	AliCentralMCCorrectionsTask::AliCentralMCCorrectionsTask(const AliCentralMCCorrectionsTask& o)
	103	: AliAnalysisTaskSE(o),
	104	fInspector(o.fInspector),
	105	fTrackDensity(),
	106	fVtxBins(0),
	107	fFirstEvent(o.fFirstEvent),
	108	fHEvents(o.fHEvents),
	109	fHEventsTr(o.fHEventsTr),
	110	fHEventsTrVtx(o.fHEventsTrVtx),
	111	fVtxAxis(10,-10,10),
	112	fEtaAxis(200,-4,6),
	113	fList(o.fList),
ea6f96fb	114	fNPhiBins(o.fNPhiBins),
ea6f96fb	115	fEffectiveCorr(o.fEffectiveCorr)
600b313f	116	{
	117	//
	118	// Copy constructor
	119	//
	120	// Parameters:
	121	// o Object to copy from
	122	//
5ca83fee	123	DGUARD(fDebug, 3,"Copy CTOR of AliCentralMCCorrectionsTask");
600b313f	124	}
	125
	126	//____________________________________________________________________
	127	AliCentralMCCorrectionsTask&
	128	AliCentralMCCorrectionsTask::operator=(const AliCentralMCCorrectionsTask& o)
	129	{
	130	//
	131	// Assignment operator
	132	//
	133	// Parameters:
	134	// o Object to assign from
	135	//
	136	// Return:
	137	// Reference to this object
	138	//
f53fb4f6	139	DGUARD(fDebug,3,"Assignment of AliCentralMCCorrectionsTask");
d015ecfe	140	if (&o == this) return *this;
600b313f	141	fInspector = o.fInspector;
	142	fTrackDensity = o.fTrackDensity;
	143	fVtxBins = o.fVtxBins;
	144	fFirstEvent = o.fFirstEvent;
	145	fHEvents = o.fHEvents;
	146	fHEventsTr = o.fHEventsTr;
	147	fHEventsTrVtx = o.fHEventsTrVtx;
	148	SetVertexAxis(o.fVtxAxis);
	149	SetEtaAxis(o.fEtaAxis);
	150	fNPhiBins = o.fNPhiBins;
ea6f96fb	151	fEffectiveCorr = o.fEffectiveCorr;
600b313f	152
	153	return *this;
	154	}
	155
	156	//____________________________________________________________________
	157	void
	158	AliCentralMCCorrectionsTask::Init()
	159	{
	160	//
	161	// Initialize the task
	162	//
	163	//
f53fb4f6	164	DGUARD(fDebug,1,"Initialize AliCentralMCCorrectionsTask");
600b313f	165	}
	166
	167	//____________________________________________________________________
	168	void
	169	AliCentralMCCorrectionsTask::SetVertexAxis(Int_t nBin, Double_t min,
	170	Double_t max)
	171	{
	172	//
	173	// Set the vertex axis to use
	174	//
	175	// Parameters:
	176	// nBins Number of bins
	177	// vzMin Least @f$z@f$ coordinate of interation point
	178	// vzMax Largest @f$z@f$ coordinate of interation point
	179	//
f53fb4f6	180	DGUARD(fDebug,3,"Set vertex axis AliCentralMCCorrectionsTask [%d,%f,%f]",
f53fb4f6	181	nBin, min, max);
bcd522ff	182	if (max < min) {
600b313f	183	Double_t tmp = min;
	184	min = max;
	185	max = tmp;
	186	}
	187	if (min < -10)
	188	AliWarning(Form("Minimum vertex %f < -10, make sure you want this",min));
	189	if (max > +10)
	190	AliWarning(Form("Minimum vertex %f > +10, make sure you want this",max));
	191	fVtxAxis.Set(nBin, min, max);
	192	}
	193	//____________________________________________________________________
	194	void
	195	AliCentralMCCorrectionsTask::SetVertexAxis(const TAxis& axis)
	196	{
	197	//
	198	// Set the vertex axis to use
	199	//
	200	// Parameters:
	201	// axis Axis
	202	//
	203	SetVertexAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax());
	204	}
	205
	206	//____________________________________________________________________
	207	void
	208	AliCentralMCCorrectionsTask::SetEtaAxis(Int_t nBin, Double_t min, Double_t max)
	209	{
	210	//
	211	// Set the eta axis to use
	212	//
	213	// Parameters:
	214	// nBins Number of bins
	215	// vzMin Least @f$\eta@f$
	216	// vzMax Largest @f$\eta@f$
	217	//
f53fb4f6	218	DGUARD(fDebug,3,"Set eta axis AliCentralMCCorrectionsTask [%d,%f,%f]",
f53fb4f6	219	nBin, min, max);
bcd522ff	220	if (max < min) {
600b313f	221	Double_t tmp = min;
	222	min = max;
	223	max = tmp;
	224	}
	225	if (min < -4)
	226	AliWarning(Form("Minimum eta %f < -4, make sure you want this",min));
	227	if (max > +6)
	228	AliWarning(Form("Minimum eta %f > +6, make sure you want this",max));
bcd522ff	229	fEtaAxis.Set(nBin, min, max);
600b313f	230	}
	231	//____________________________________________________________________
	232	void
	233	AliCentralMCCorrectionsTask::SetEtaAxis(const TAxis& axis)
	234	{
	235	//
	236	// Set the eta axis to use
	237	//
	238	// Parameters:
	239	// axis Axis
	240	//
	241	SetEtaAxis(axis.GetNbins(),axis.GetXmin(),axis.GetXmax());
	242	}
	243
	244	//____________________________________________________________________
	245	void
	246	AliCentralMCCorrectionsTask::DefineBins(TList* l)
	247	{
f53fb4f6	248	DGUARD(fDebug,1,"Define bins in AliCentralMCCorrectionsTask");
600b313f	249	if (!fVtxBins) fVtxBins = new TObjArray(fVtxAxis.GetNbins(), 1);
	250	if (fVtxBins->GetEntries() > 0) return;
	251
	252	fVtxBins->SetOwner();
	253	for (Int_t i = 1; i <= fVtxAxis.GetNbins(); i++) {
	254	Double_t low = fVtxAxis.GetBinLowEdge(i);
	255	Double_t high = fVtxAxis.GetBinUpEdge(i);
	256	VtxBin* bin = new VtxBin(low, high, fEtaAxis, fNPhiBins);
	257	fVtxBins->AddAt(bin, i);
5934a3e3	258	bin->CreateOutputObjects(l);
600b313f	259	}
	260	}
	261
	262	//____________________________________________________________________
	263	void
	264	AliCentralMCCorrectionsTask::UserCreateOutputObjects()
	265	{
	266	//
	267	// Create output objects
	268	//
	269	//
f53fb4f6	270	DGUARD(fDebug,1,"Create user output for AliCentralMCCorrectionsTask");
600b313f	271	fList = new TList;
	272	fList->SetOwner();
	273	fList->SetName(Form("%sSums", GetName()));
	274
	275	DefineBins(fList);
	276
	277	fHEvents = new TH1I(GetEventName(false,false),
	278	"Number of all events",
	279	fVtxAxis.GetNbins(),
	280	fVtxAxis.GetXmin(),
	281	fVtxAxis.GetXmax());
	282	fHEvents->SetXTitle("v_{z} [cm]");
	283	fHEvents->SetYTitle("# of events");
	284	fHEvents->SetFillColor(kBlue+1);
	285	fHEvents->SetFillStyle(3001);
	286	fHEvents->SetDirectory(0);
	287	fList->Add(fHEvents);
	288
	289	fHEventsTr = new TH1I(GetEventName(true, false),
	290	"Number of triggered events",
	291	fVtxAxis.GetNbins(),
	292	fVtxAxis.GetXmin(),
	293	fVtxAxis.GetXmax());
	294	fHEventsTr->SetXTitle("v_{z} [cm]");
	295	fHEventsTr->SetYTitle("# of events");
	296	fHEventsTr->SetFillColor(kRed+1);
	297	fHEventsTr->SetFillStyle(3001);
	298	fHEventsTr->SetDirectory(0);
	299	fList->Add(fHEventsTr);
	300
	301	fHEventsTrVtx = new TH1I(GetEventName(true,true),
	302	"Number of events w/trigger and vertex",
	303	fVtxAxis.GetNbins(),
	304	fVtxAxis.GetXmin(),
	305	fVtxAxis.GetXmax());
	306	fHEventsTrVtx->SetXTitle("v_{z} [cm]");
	307	fHEventsTrVtx->SetYTitle("# of events");
	308	fHEventsTrVtx->SetFillColor(kBlue+1);
	309	fHEventsTrVtx->SetFillStyle(3001);
	310	fHEventsTrVtx->SetDirectory(0);
	311	fList->Add(fHEventsTrVtx);
	312
	313	// Copy axis objects to output
	314	TH1* vtxAxis = new TH1D("vtxAxis", "Vertex axis",
	315	fVtxAxis.GetNbins(),
	316	fVtxAxis.GetXmin(),
	317	fVtxAxis.GetXmax());
	318	TH1* etaAxis = new TH1D("etaAxis", "Eta axis",
	319	fEtaAxis.GetNbins(),
	320	fEtaAxis.GetXmin(),
	321	fEtaAxis.GetXmax());
	322	fList->Add(vtxAxis);
	323	fList->Add(etaAxis);
	324
	325	AliInfo(Form("Initialising sub-routines: %p, %p",
	326	&fInspector, &fTrackDensity));
5934a3e3	327	fInspector.CreateOutputObjects(fList);
5934a3e3	328	fTrackDensity.CreateOutputObjects(fList);
600b313f	329
	330	PostData(1, fList);
	331	}
	332	//____________________________________________________________________
	333	void
	334	AliCentralMCCorrectionsTask::UserExec(Option_t*)
	335	{
	336	//
	337	// Process each event
	338	//
	339	// Parameters:
	340	// option Not used
	341	//
	342
f53fb4f6	343	DGUARD(fDebug,1,"AliCentralMCCorrectionsTask process an event");
600b313f	344	// Get the input data - MC event
	345	AliMCEvent* mcEvent = MCEvent();
	346	if (!mcEvent) {
	347	AliWarning("No MC event found");
	348	return;
	349	}
	350
	351	// Get the input data - ESD event
	352	AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
	353	if (!esd) {
	354	AliWarning("No ESD event found for input event");
	355	return;
	356	}
	357
	358	//--- Read run information -----------------------------------------
	359	if (fFirstEvent && esd->GetESDRun()) {
	360	fInspector.ReadRunDetails(esd);
	361
	362	AliInfo(Form("Initializing with parameters from the ESD:\n"
	363	" AliESDEvent::GetBeamEnergy() ->%f\n"
	364	" AliESDEvent::GetBeamType() ->%s\n"
	365	" AliESDEvent::GetCurrentL3() ->%f\n"
	366	" AliESDEvent::GetMagneticField()->%f\n"
	367	" AliESDEvent::GetRunNumber() ->%d\n",
	368	esd->GetBeamEnergy(),
	369	esd->GetBeamType(),
	370	esd->GetCurrentL3(),
	371	esd->GetMagneticField(),
	372	esd->GetRunNumber()));
	373
5934a3e3	374	fInspector.SetupForData(fVtxAxis);
576472c1	375
600b313f	376	Print();
	377	fFirstEvent = false;
	378	}
	379
	380	// Some variables
3e425989	381	UInt_t triggers = 0; // Trigger bits
	382	Bool_t lowFlux = true; // Low flux flag
	383	UShort_t iVz = 0; // Vertex bin from ESD
5ca83fee	384	TVector3 ip; // Z coordinate from ESD
3e425989	385	Double_t cent = -1; // Centrality
	386	UShort_t iVzMc = 0; // Vertex bin from MC
	387	Double_t vZMc = 1000; // Z coordinate of IP vertex from MC
	388	Double_t b = -1; // Impact parameter
	389	Double_t cMC = -1; // Centrality estimate from b
	390	Int_t nPart = -1; // Number of participants
	391	Int_t nBin = -1; // Number of binary collisions
	392	Double_t phiR = 1000; // Reaction plane from MC
	393	UShort_t nClusters = 0; // Number of SPD clusters
600b313f	394	// Process the data
5ca83fee	395	UInt_t retESD = fInspector.Process(esd, triggers, lowFlux, iVz, ip,
600b313f	396	cent, nClusters);
600b313f	397	fInspector.ProcessMC(mcEvent, triggers, iVzMc, vZMc,
241cca4d	398	b, cMC, nPart, nBin, phiR);
600b313f	399
	400	Bool_t isInel = triggers & AliAODForwardMult::kInel;
	401	Bool_t hasVtx = retESD == AliFMDMCEventInspector::kOk;
	402
	403	// Fill the event count histograms
	404	if (isInel) fHEventsTr->Fill(vZMc);
	405	if (isInel && hasVtx) fHEventsTrVtx->Fill(vZMc);
	406	fHEvents->Fill(vZMc);
	407
	408	// Now find our vertex bin object
	409	VtxBin* bin = 0;
	410	if (iVzMc > 0 && iVzMc <= fVtxAxis.GetNbins())
	411	bin = static_cast<VtxBin*>(fVtxBins->At(iVzMc));
	412	if (!bin) {
	413	// AliError(Form("No vertex bin object @ %d (%f)", iVzMc, vZMc));
	414	return;
	415	}
	416
	417	// Now process our input data and store in own ESD object
	418	fTrackDensity.Calculate(mcEvent, vZMc, bin->fHits, bin->fPrimary);
ea6f96fb	419
	420	// Get the ESD object
	421	const AliMultiplicity* spdmult = esd->GetMultiplicity();
	422
	423	// Count number of tracklets per bin
	424	for(Int_t j = 0; j< spdmult->GetNumberOfTracklets();j++)
	425	bin->fClusters->Fill(spdmult->GetEta(j),spdmult->GetPhi(j));
	426	//...and then the unused clusters in layer 1
	427	for(Int_t j = 0; j< spdmult->GetNumberOfSingleClusters();j++) {
	428	Double_t eta = -TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.));
	429	bin->fClusters->Fill(eta, spdmult->GetPhiSingle(j));
	430	}
	431
	432	// Count events
600b313f	433	bin->fCounts->Fill(0.5);
bcd522ff	434
bcd522ff	435	PostData(1, fList);
600b313f	436	}
	437
	438	//____________________________________________________________________
	439	void
	440	AliCentralMCCorrectionsTask::Terminate(Option_t*)
	441	{
	442	//
	443	// End of job
	444	//
	445	// Parameters:
	446	// option Not used
	447	//
f53fb4f6	448	DGUARD(fDebug,1,"AliCentralMCCorrectionsTask analyse merged output");
600b313f	449
	450	fList = dynamic_cast<TList*>(GetOutputData(1));
	451	if (!fList) {
	452	AliError("No output list defined");
	453	return;
	454	}
	455
bcd522ff	456	DefineBins(fList);
bcd522ff	457
600b313f	458	// Output list
	459	TList* output = new TList;
	460	output->SetOwner();
	461	output->SetName(Form("%sResults", GetName()));
	462
	463	// --- Fill correction object --------------------------------------
	464	AliCentralCorrSecondaryMap* corr = new AliCentralCorrSecondaryMap;
	465	corr->SetVertexAxis(fVtxAxis);
	466	// corr->SetEtaAxis(fEtaAxis);
	467
438d5c4e	468	AliCentralCorrAcceptance* acorr = new AliCentralCorrAcceptance;
ea6f96fb	469	acorr->SetVertexAxis(fVtxAxis);
ea6f96fb	470
600b313f	471	TIter next(fVtxBins);
	472	VtxBin* bin = 0;
	473	UShort_t iVz = 1;
	474	while ((bin = static_cast<VtxBin*>(next())))
5934a3e3	475	bin->Terminate(fList, output, iVz++, fEffectiveCorr, corr,acorr);
600b313f	476
600b313f	477	output->Add(corr);
438d5c4e	478	output->Add(acorr);
600b313f	479
	480	PostData(2, output);
	481	}
	482
	483	//____________________________________________________________________
	484	void
	485	AliCentralMCCorrectionsTask::Print(Option_t* option) const
	486	{
	487	std::cout << ClassName() << "\n"
	488	<< " Vertex bins: " << fVtxAxis.GetNbins() << '\n'
	489	<< " Vertex range: [" << fVtxAxis.GetXmin()
	490	<< "," << fVtxAxis.GetXmax() << "]\n"
	491	<< " Eta bins: " << fEtaAxis.GetNbins() << '\n'
	492	<< " Eta range: [" << fEtaAxis.GetXmin()
	493	<< "," << fEtaAxis.GetXmax() << "]\n"
	494	<< " # of phi bins: " << fNPhiBins
	495	<< std::endl;
	496	gROOT->IncreaseDirLevel();
	497	fInspector.Print(option);
	498	fTrackDensity.Print(option);
	499	gROOT->DecreaseDirLevel();
	500	}
	501
	502	//====================================================================
	503	const char*
	504	AliCentralMCCorrectionsTask::VtxBin::BinName(Double_t low,
	505	Double_t high)
	506	{
	507	static TString buf;
	508	buf = Form("vtx%+05.1f_%+05.1f", low, high);
	509	buf.ReplaceAll("+", "p");
	510	buf.ReplaceAll("-", "m");
	511	buf.ReplaceAll(".", "d");
	512	return buf.Data();
	513	}
	514
	515
	516	//____________________________________________________________________
	517	AliCentralMCCorrectionsTask::VtxBin::VtxBin()
	518	: fHits(0),
ea6f96fb	519	fClusters(0),
600b313f	520	fPrimary(0),
	521	fCounts(0)
	522	{
	523	}
	524	//____________________________________________________________________
	525	AliCentralMCCorrectionsTask::VtxBin::VtxBin(Double_t low,
	526	Double_t high,
	527	const TAxis& axis,
	528	UShort_t nPhi)
	529	: TNamed(BinName(low, high),
	530	Form("%+5.1fcm<v_{z}<%+5.1fcm", low, high)),
	531	fHits(0),
ea6f96fb	532	fClusters(0),
600b313f	533	fPrimary(0),
	534	fCounts(0)
	535	{
	536	fPrimary = new TH2D("primary", "Primaries",
	537	axis.GetNbins(), axis.GetXmin(), axis.GetXmax(),
	538	nPhi, 0, 2*TMath::Pi());
	539	fPrimary->SetXTitle("#eta");
	540	fPrimary->SetYTitle("#varphi [radians]");
	541	fPrimary->Sumw2();
	542	fPrimary->SetDirectory(0);
	543
	544	fHits = static_cast<TH2D*>(fPrimary->Clone("hits"));
	545	fHits->SetTitle("Hits");
	546	fHits->SetDirectory(0);
	547
ea6f96fb	548	fClusters = static_cast<TH2D*>(fPrimary->Clone("clusters"));
	549	fClusters->SetTitle("Clusters");
	550	fClusters->SetDirectory(0);
	551
600b313f	552	fCounts = new TH1D("counts", "Counts", 1, 0, 1);
	553	fCounts->SetXTitle("Events");
	554	fCounts->SetYTitle("# of Events");
	555	fCounts->SetDirectory(0);
	556	}
	557
	558	//____________________________________________________________________
	559	AliCentralMCCorrectionsTask::VtxBin::VtxBin(const VtxBin& o)
	560	: TNamed(o),
	561	fHits(0),
ea6f96fb	562	fClusters(0),
600b313f	563	fPrimary(0),
	564	fCounts(0)
	565	{
	566	if (o.fHits) {
	567	fHits = static_cast<TH2D*>(o.fHits->Clone());
	568	fHits->SetDirectory(0);
	569	}
ea6f96fb	570	if (o.fClusters) {
	571	fClusters = static_cast<TH2D*>(o.fClusters->Clone());
	572	fClusters->SetDirectory(0);
	573	}
600b313f	574	if (o.fPrimary) {
	575	fPrimary = static_cast<TH2D*>(o.fPrimary->Clone());
	576	fPrimary->SetDirectory(0);
	577	}
	578	if (o.fCounts) {
	579	fCounts = static_cast<TH1D*>(o.fCounts->Clone());
	580	fCounts->SetDirectory(0);
	581	}
	582	}
	583
	584	//____________________________________________________________________
	585	AliCentralMCCorrectionsTask::VtxBin&
	586	AliCentralMCCorrectionsTask::VtxBin::operator=(const VtxBin& o)
	587	{
d015ecfe	588	if (&o == this) return *this;
600b313f	589	TNamed::operator=(o);
ea6f96fb	590	fHits = 0;
	591	fPrimary = 0;
	592	fClusters = 0;
	593	fCounts = 0;
600b313f	594	if (o.fHits) {
	595	fHits = static_cast<TH2D*>(o.fHits->Clone());
	596	fHits->SetDirectory(0);
	597	}
ea6f96fb	598	if (o.fClusters) {
	599	fClusters = static_cast<TH2D*>(o.fClusters->Clone());
	600	fClusters->SetDirectory(0);
	601	}
600b313f	602	if (o.fPrimary) {
	603	fPrimary = static_cast<TH2D*>(o.fPrimary->Clone());
	604	fPrimary->SetDirectory(0);
	605	}
	606	if (o.fCounts) {
	607	fCounts = static_cast<TH1D*>(o.fCounts->Clone());
	608	fCounts->SetDirectory(0);
	609	}
	610	return *this;
	611	}
	612
	613	//____________________________________________________________________
	614	void
5934a3e3	615	AliCentralMCCorrectionsTask::VtxBin::CreateOutputObjects(TList* l)
600b313f	616	{
	617	TList* d = new TList;
	618	d->SetName(GetName());
	619	d->SetOwner();
	620	l->Add(d);
	621
	622	d->Add(fHits);
ea6f96fb	623	d->Add(fClusters);
600b313f	624	d->Add(fPrimary);
	625	d->Add(fCounts);
	626	}
	627	//____________________________________________________________________
	628	void
5934a3e3	629	AliCentralMCCorrectionsTask::VtxBin::Terminate(const TList* input,
600b313f	630	TList* output,
600b313f	631	UShort_t iVz,
ea6f96fb	632	Bool_t effectiveCorr,
	633	AliCentralCorrSecondaryMap* map,
	634	AliCentralCorrAcceptance* acorr)
600b313f	635	{
	636	TList* out = new TList;
	637	out->SetName(GetName());
	638	out->SetOwner();
	639	output->Add(out);
	640
	641	TList* l = static_cast<TList*>(input->FindObject(GetName()));
	642	if (!l) {
	643	AliError(Form("List %s not found in %s", GetName(), input->GetName()));
	644	return;
	645	}
	646
	647
	648	TH2D* hits = static_cast<TH2D*>(l->FindObject("hits"));
ea6f96fb	649	TH2D* clus = static_cast<TH2D*>(l->FindObject("clusters"));
600b313f	650	TH2D* prim = static_cast<TH2D*>(l->FindObject("primary"));
	651	if (!hits \|\| !prim) {
	652	AliError(Form("Missing histograms: %p, %p", hits, prim));
	653	return;
	654	}
	655
ea6f96fb	656	TH2D* h = 0;
	657	if (effectiveCorr) h = static_cast<TH2D*>(clus->Clone("bgCorr"));
	658	else h = static_cast<TH2D*>(hits->Clone("bgCorr"));
600b313f	659	h->SetDirectory(0);
984a7263	660	h->Divide(prim);
600b313f	661
ea6f96fb	662	TH1D* acc = new TH1D(Form("SPDacc_vrtbin_%d",iVz),
438d5c4e	663	"Acceptance correction for SPD" ,
	664	fPrimary->GetXaxis()->GetNbins(),
	665	fPrimary->GetXaxis()->GetXmin(),
	666	fPrimary->GetXaxis()->GetXmax());
	667	TH1F* accden = static_cast<TH1F*>(acc->Clone(Form("%s_den",acc->GetName())));
ea6f96fb	668
	669	for(Int_t xx = 1; xx <=h->GetNbinsX(); xx++) {
	670	for(Int_t yy = 1; yy <=h->GetNbinsY(); yy++) {
	671	if(TMath::Abs(h->GetXaxis()->GetBinCenter(xx)) > 1.9) {
	672	h->SetBinContent(xx,yy,0.);
	673	h->SetBinError(xx,yy,0.);
	674	}
	675	if(h->GetBinContent(xx,yy) > 0.9)
	676	acc->Fill(h->GetXaxis()->GetBinCenter(xx));
	677	else {
	678	h->SetBinContent(xx,yy,0.);
	679	h->SetBinError(xx,yy,0.);
	680	}
	681	accden->Fill(h->GetXaxis()->GetBinCenter(xx));
ea6f96fb	682	}
	683	}
	684	acc->Divide(accden);
600b313f	685
ea6f96fb	686	map->SetCorrection(iVz, h);
ea6f96fb	687	acorr->SetCorrection(iVz, acc);
600b313f	688	out->Add(h);
ea6f96fb	689	out->Add(acc);
600b313f	690	}
	691
	692	//
	693	// EOF
	694	//