[u/mrichter/AliRoot.git] / PWG2 / 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 <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" : ""));
  }
  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(),
    fTrackDensity(),
    fVtxBins(0),
    fFirstEvent(true),
    fHEvents(0), 
    fHEventsTr(0), 
    fHEventsTrVtx(0),
    fVtxAxis(),
    fEtaAxis(),
    fList(),
    fNPhiBins(40)    
{
  // 
  // Constructor 
  // 
  // Parameters:
  //    name Name of task 
  //
}

//____________________________________________________________________
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(40)    
{
  // 
  // Constructor 
  // 
  // 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 
  //
  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 
  // 
  //
}

//____________________________________________________________________
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) max = -min;
  if (min < max) { 
    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$ 
  //
  if (max < min) max = -min;
  if (min < max) { 
    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));
  fVtxAxis.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)
{
  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 
  // 
  //
  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
  //  

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

//____________________________________________________________________
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;
  }

  // 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)
{
  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(fPrimary);

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