[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 "AliForwardUtil.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)    
{
  // 
  // Constructor 
  // 
  // Parameters:
  //    name Name of task 
  //
  DGUARD(fDebug,0,"Default construction 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)    
{
  // 
  // Constructor 
  // 
  // Parameters:
  //    name Name of task 
  //
  DGUARD(fDebug,0,"Named construction of AliCentralMCCorrectionsTask: %s",name);
  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 
  //
  DGUARD(fDebug,0,"Copy construction 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;

  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->DefineOutput(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.DefineOutput(fList);
  fInspector.Init(fVtxAxis);
  fTrackDensity.DefineOutput(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()));

    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
  Double_t cMC;      // Centrality estimate from b
  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, 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);
  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);

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

//____________________________________________________________________
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), 
    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), 
    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);

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

//____________________________________________________________________
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* d = new TList;
  d->SetName(GetName());
  d->SetOwner();
  l->Add(d);

  d->Add(fHits);
  d->Add(fPrimary);
  d->Add(fCounts);
}
//____________________________________________________________________
void
AliCentralMCCorrectionsTask::VtxBin::Finish(const TList* input, 
					    TList* output, 
					    UShort_t iVz, 
					    AliCentralCorrSecondaryMap* map)
{
  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*   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);

  out->Add(h);
}

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