[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliCentralMultiplicityTask.cxx

//====================================================================
// 
// Base class for classes that calculate the multiplicity in the
// central region event-by-event
// 
// Inputs: 
//   - AliESDEvent 
//
// Outputs: 
//   - AliAODCentralMult 
// 
// Histograms 
//   
// Corrections used 
#include "AliCentralMultiplicityTask.h"
#include "AliForwardCorrectionManager.h"
#include "AliForwardUtil.h"
#include "AliLog.h"
#include "AliAODHandler.h"
#include "AliInputEventHandler.h"
#include "AliESDInputHandler.h"
#include "AliAnalysisManager.h"
#include "AliESDEvent.h"
#include "AliMultiplicity.h"
#include <TROOT.h>
#include <TFile.h>
#include <TError.h>
#include <iostream>
#include <iomanip>

//====================================================================
AliCentralMultiplicityTask::AliCentralMultiplicityTask(const char* name) 
  : AliAnalysisTaskSE(name),
    fData(0),
    fList(0),
    fAODCentral(kFALSE),
    fManager()
{
  
  DefineOutput(1, TList::Class());
}
//____________________________________________________________________
void AliCentralMultiplicityTask::UserCreateOutputObjects() 
{

  AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
  AliAODHandler*      ah = 
    dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
  if (!ah) AliFatal("No AOD output handler set in analysis manager");
  
  
  TObject* obj = &fAODCentral;
  ah->AddBranch("AliAODCentralMult", &obj);
  
  fList = new TList();
  PostData(1,fList);
  
}
//____________________________________________________________________
void AliCentralMultiplicityTask::UserExec(Option_t* /*option*/) 
{
  
  AliESDInputHandler* eventHandler = 
    dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()
				       ->GetInputEventHandler());
  if (!eventHandler) {
    AliWarning("No inputhandler found for this event!");
    return;}
  
  AliESDEvent* esd = eventHandler->GetEvent();
  
  //Selecting only events with |valid vertex| < 10 cm
  const AliESDVertex* vertex = esd->GetPrimaryVertexSPD();
  if (!vertex) return;
  if(!(vertex->GetStatus())) return;
  if(vertex->GetNContributors() <= 0) return ;
  if(vertex->GetZRes() > 0.1 ) return;
  Double_t vertexXYZ[3]={0,0,0};
  vertex->GetXYZ(vertexXYZ);
  if(TMath::Abs(vertexXYZ[2]) > 10) return;
  
  Double_t delta           = 2 ;
  Double_t vertexBinDouble = (vertexXYZ[2] + 10) / delta;
  //HHD: The vtxbins are 1-10, not 0-9
  Int_t    vtxbin          = Int_t(vertexBinDouble + 1) ; 
  
  // Make sure AOD is filled
  AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
  AliAODHandler*      ah = 
    dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
  if (!ah)  
    AliFatal("No AOD output handler set in analysis manager");
  
  ah->SetFillAOD(kTRUE);
  
  //Doing analysis
  fAODCentral.Clear("");
  TH2D *aodHist = &(fAODCentral.GetHistogram());
  
  const AliMultiplicity* spdmult = esd->GetMultiplicity();
  //Filling clusters in layer 1 used for tracklets...
  for(Int_t j = 0; j< spdmult->GetNumberOfTracklets();j++)
    aodHist->Fill(spdmult->GetEta(j),spdmult->GetPhi(j));

  //...and then the unused ones in layer 1 
  for(Int_t j = 0; j< spdmult->GetNumberOfSingleClusters();j++) 
    aodHist->Fill(-TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.)),
		  spdmult->GetPhiSingle(j));
  
  
  // Corrections
  TH2D* hSecMap     = fManager.GetSecMapCorrection(vtxbin);
  TH1D* hAcceptance = fManager.GetAcceptanceCorrection(vtxbin);
  
  aodHist->Divide(hSecMap);
  
  for(Int_t nx = 1; nx <= aodHist->GetNbinsX(); nx++) {
    Float_t acccor = hAcceptance->GetBinContent(nx);
    
    Bool_t etabinSeen = kFALSE;  
    for(Int_t ny = 1; ny <= aodHist->GetNbinsY(); ny++) {
      Float_t aodvalue = aodHist->GetBinContent(nx,ny);
      Float_t seccor = hSecMap->GetBinContent(nx,ny);
      if(seccor > 0.5) etabinSeen = kTRUE;
      if(aodvalue < 0.000001) { aodHist->SetBinContent(nx,ny, 0); continue; }
      
      Float_t aodnew   = aodvalue / acccor ;
      aodHist->SetBinContent(nx,ny, aodnew);
      Float_t aodErr   = aodHist->GetBinError(nx,ny);
      Float_t accErr   = hAcceptance->GetBinError(nx);
      Float_t error    = aodnew *TMath::Sqrt(TMath::Power(aodErr/aodvalue,2) +
					     TMath::Power(accErr/acccor,2) );
      aodHist->SetBinError(nx,ny,error);
      
    }
    //Filling underflow bin if we eta bin is in range
    if(etabinSeen) aodHist->SetBinContent(nx,0, 1.);
  }  

  PostData(1,fList);
}
//____________________________________________________________________
void AliCentralMultiplicityTask::Terminate(Option_t* /*option*/) 
{
}
//____________________________________________________________________
void
AliCentralMultiplicityTask::Print(Option_t* /*option*/) const
{
}
//====================================================================
AliCentralMultiplicityTask::Manager::Manager() :
  fAcceptancePath("$ALICE_ROOT/PWG2/FORWARD/corrections/CentralAcceptance"),
  fSecMapPath("$ALICE_ROOT/PWG2/FORWARD/corrections/CentralSecMap"),
  fAcceptance(),
  fSecmap(),
  fAcceptanceName("centralacceptance"),
  fSecMapName("centralsecmap")
{


}
//____________________________________________________________________
AliCentralMultiplicityTask::Manager::Manager(const Manager& o) 
  :fAcceptancePath(o.fAcceptancePath),
   fSecMapPath(o.fSecMapPath),
   fAcceptance(o.fAcceptance),
   fSecmap(o.fSecmap),
   fAcceptanceName(o.fAcceptanceName),
   fSecMapName(o.fSecMapName) 
{}
//____________________________________________________________________
AliCentralMultiplicityTask::Manager&
AliCentralMultiplicityTask::Manager::operator=(const Manager& o)
{
  fAcceptancePath = o.fAcceptancePath;
  fSecMapPath     = o.fSecMapPath;
  fAcceptance     = o.fAcceptance;
  fSecmap         = o.fSecmap;
  fAcceptanceName = o.fAcceptanceName;
  fSecMapName     = o.fSecMapName;
  return *this;
}

//____________________________________________________________________
const char* 
AliCentralMultiplicityTask::Manager::GetFullFileName(UShort_t what, 
						     UShort_t sys, 
						     UShort_t sNN,  
						     Short_t  field) const
{
  return Form("%s/%s",
	      what == 0 ? GetSecMapPath() : GetAcceptancePath(), 
	      GetFileName(what, sys, sNN, field));
}

//____________________________________________________________________
const char* 
AliCentralMultiplicityTask::Manager::GetFileName(UShort_t  what ,
						 UShort_t  sys, 
						 UShort_t  sNN,
						 Short_t   field) const
{
  // Must be static - otherwise the data may disappear on return from
  // this member function
  static TString fname = "";
  
  switch(what) {
  case 0:  fname.Append(fSecMapName.Data());     break;
  case 1:  fname.Append(fAcceptanceName.Data()); break;
  default:
    ::Error("GetFileName", 
	    "Invalid indentifier %d for central object, must be 0 or 1!", what);
    break;
  }
  fname.Append(Form("_%s_%04dGeV_%c%1dkG.root", 
		    AliForwardUtil::CollisionSystemString(sys), 
		    sNN, (field < 0 ? 'm' : 'p'), TMath::Abs(field)));
  
  return fname.Data();
}

//____________________________________________________________________
TH2D* 
AliCentralMultiplicityTask::Manager::GetSecMapCorrection(UShort_t vtxbin) const
{
  if (!fSecmap) { 
    ::Warning("GetSecMapCorrection","No secondary map defined");
    return 0;
  }
  return fSecmap->GetCorrection(vtxbin);
}
//____________________________________________________________________
TH1D* 
AliCentralMultiplicityTask::Manager::GetAcceptanceCorrection(UShort_t vtxbin) 
  const 
{
  if (!fAcceptance) { 
    ::Warning("GetAcceptanceCorrection","No acceptance map defined");
    return 0;
  }
  return fAcceptance->GetCorrection(vtxbin);
}

//____________________________________________________________________
void 
AliCentralMultiplicityTask::Manager::Init(UShort_t  sys, 
					  UShort_t  sNN,
					  Short_t   field) 
{
  TFile fsec(GetFullFileName(0,sys,sNN,field));
  fSecmap = 
    dynamic_cast<AliCentralCorrSecondaryMap*>(fsec.Get(fSecMapName.Data()));  
  if(!fSecmap) {
    ::Error("Init", "no central Secondary Map found!") ;
    return;
  }
  TFile facc(GetFullFileName(1,sys,sNN,field));
  fAcceptance = 
    dynamic_cast<AliCentralCorrAcceptance*>(facc.Get(fAcceptanceName.Data()));
  if(!fAcceptance) {
    ::Error("Init", "no central Acceptance found!") ;
    return;
  }

}
//
// EOF
//
Commit	Line	Data
6f791cc3	1	//====================================================================
	2	//
	3	// Base class for classes that calculate the multiplicity in the
	4	// central region event-by-event
	5	//
	6	// Inputs:
	7	// - AliESDEvent
	8	//
	9	// Outputs:
	10	// - AliAODCentralMult
	11	//
	12	// Histograms
	13	//
	14	// Corrections used
	15	#include "AliCentralMultiplicityTask.h"
	16	#include "AliForwardCorrectionManager.h"
	17	#include "AliForwardUtil.h"
	18	#include "AliLog.h"
	19	#include "AliAODHandler.h"
	20	#include "AliInputEventHandler.h"
	21	#include "AliESDInputHandler.h"
	22	#include "AliAnalysisManager.h"
	23	#include "AliESDEvent.h"
	24	#include "AliMultiplicity.h"
	25	#include <TROOT.h>
	26	#include <TFile.h>
3e478dba	27	#include <TError.h>
6f791cc3	28	#include <iostream>
	29	#include <iomanip>
	30
	31	//====================================================================
	32	AliCentralMultiplicityTask::AliCentralMultiplicityTask(const char* name)
	33	: AliAnalysisTaskSE(name),
	34	fData(0),
	35	fList(0),
	36	fAODCentral(kFALSE),
	37	fManager()
	38	{
	39
	40	DefineOutput(1, TList::Class());
	41	}
	42	//____________________________________________________________________
3e478dba	43	void AliCentralMultiplicityTask::UserCreateOutputObjects()
3e478dba	44	{
6f791cc3	45
	46	AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
	47	AliAODHandler* ah =
	48	dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
	49	if (!ah) AliFatal("No AOD output handler set in analysis manager");
	50
	51
	52	TObject* obj = &fAODCentral;
	53	ah->AddBranch("AliAODCentralMult", &obj);
	54
	55	fList = new TList();
	56	PostData(1,fList);
	57
	58	}
	59	//____________________________________________________________________
3e478dba	60	void AliCentralMultiplicityTask::UserExec(Option_t* /option/)
3e478dba	61	{
6f791cc3	62
	63	AliESDInputHandler* eventHandler =
	64	dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()
	65	->GetInputEventHandler());
	66	if (!eventHandler) {
	67	AliWarning("No inputhandler found for this event!");
	68	return;}
	69
	70	AliESDEvent* esd = eventHandler->GetEvent();
	71
	72	//Selecting only events with \|valid vertex\| < 10 cm
	73	const AliESDVertex* vertex = esd->GetPrimaryVertexSPD();
	74	if (!vertex) return;
	75	if(!(vertex->GetStatus())) return;
	76	if(vertex->GetNContributors() <= 0) return ;
	77	if(vertex->GetZRes() > 0.1 ) return;
	78	Double_t vertexXYZ[3]={0,0,0};
	79	vertex->GetXYZ(vertexXYZ);
	80	if(TMath::Abs(vertexXYZ[2]) > 10) return;
	81
3e478dba	82	Double_t delta = 2 ;
6f791cc3	83	Double_t vertexBinDouble = (vertexXYZ[2] + 10) / delta;
3e478dba	84	//HHD: The vtxbins are 1-10, not 0-9
3e478dba	85	Int_t vtxbin = Int_t(vertexBinDouble + 1) ;
6f791cc3	86
	87	// Make sure AOD is filled
	88	AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
	89	AliAODHandler* ah =
	90	dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
	91	if (!ah)
	92	AliFatal("No AOD output handler set in analysis manager");
	93
	94	ah->SetFillAOD(kTRUE);
	95
	96	//Doing analysis
	97	fAODCentral.Clear("");
	98	TH2D *aodHist = &(fAODCentral.GetHistogram());
	99
	100	const AliMultiplicity* spdmult = esd->GetMultiplicity();
	101	//Filling clusters in layer 1 used for tracklets...
3e478dba	102	for(Int_t j = 0; j< spdmult->GetNumberOfTracklets();j++)
6f791cc3	103	aodHist->Fill(spdmult->GetEta(j),spdmult->GetPhi(j));
3e478dba	104
6f791cc3	105	//...and then the unused ones in layer 1
3e478dba	106	for(Int_t j = 0; j< spdmult->GetNumberOfSingleClusters();j++)
6f791cc3	107	aodHist->Fill(-TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.)),
3e478dba	108	spdmult->GetPhiSingle(j));
6f791cc3	109
6f791cc3	110
3e478dba	111	// Corrections
6f791cc3	112	TH2D* hSecMap = fManager.GetSecMapCorrection(vtxbin);
	113	TH1D* hAcceptance = fManager.GetAcceptanceCorrection(vtxbin);
	114
	115	aodHist->Divide(hSecMap);
	116
	117	for(Int_t nx = 1; nx <= aodHist->GetNbinsX(); nx++) {
	118	Float_t acccor = hAcceptance->GetBinContent(nx);
	119
	120	Bool_t etabinSeen = kFALSE;
	121	for(Int_t ny = 1; ny <= aodHist->GetNbinsY(); ny++) {
	122	Float_t aodvalue = aodHist->GetBinContent(nx,ny);
	123	Float_t seccor = hSecMap->GetBinContent(nx,ny);
	124	if(seccor > 0.5) etabinSeen = kTRUE;
	125	if(aodvalue < 0.000001) { aodHist->SetBinContent(nx,ny, 0); continue; }
	126
	127	Float_t aodnew = aodvalue / acccor ;
	128	aodHist->SetBinContent(nx,ny, aodnew);
3e478dba	129	Float_t aodErr = aodHist->GetBinError(nx,ny);
	130	Float_t accErr = hAcceptance->GetBinError(nx);
	131	Float_t error = aodnew *TMath::Sqrt(TMath::Power(aodErr/aodvalue,2) +
	132	TMath::Power(accErr/acccor,2) );
6f791cc3	133	aodHist->SetBinError(nx,ny,error);
	134
	135	}
	136	//Filling underflow bin if we eta bin is in range
	137	if(etabinSeen) aodHist->SetBinContent(nx,0, 1.);
6f791cc3	138	}
	139
	140	PostData(1,fList);
6f791cc3	141	}
6f791cc3	142	//____________________________________________________________________
3e478dba	143	void AliCentralMultiplicityTask::Terminate(Option_t* /option/)
3e478dba	144	{
6f791cc3	145	}
	146	//____________________________________________________________________
	147	void
	148	AliCentralMultiplicityTask::Print(Option_t* /option/) const
	149	{
6f791cc3	150	}
3e478dba	151	//====================================================================
6f791cc3	152	AliCentralMultiplicityTask::Manager::Manager() :
	153	fAcceptancePath("$ALICE_ROOT/PWG2/FORWARD/corrections/CentralAcceptance"),
	154	fSecMapPath("$ALICE_ROOT/PWG2/FORWARD/corrections/CentralSecMap"),
	155	fAcceptance(),
	156	fSecmap(),
	157	fAcceptanceName("centralacceptance"),
	158	fSecMapName("centralsecmap")
	159	{
	160
	161
	162	}
	163	//____________________________________________________________________
3e478dba	164	AliCentralMultiplicityTask::Manager::Manager(const Manager& o)
	165	:fAcceptancePath(o.fAcceptancePath),
	166	fSecMapPath(o.fSecMapPath),
	167	fAcceptance(o.fAcceptance),
	168	fSecmap(o.fSecmap),
	169	fAcceptanceName(o.fAcceptanceName),
	170	fSecMapName(o.fSecMapName)
	171	{}
	172	//____________________________________________________________________
	173	AliCentralMultiplicityTask::Manager&
	174	AliCentralMultiplicityTask::Manager::operator=(const Manager& o)
	175	{
	176	fAcceptancePath = o.fAcceptancePath;
	177	fSecMapPath = o.fSecMapPath;
	178	fAcceptance = o.fAcceptance;
	179	fSecmap = o.fSecmap;
	180	fAcceptanceName = o.fAcceptanceName;
	181	fSecMapName = o.fSecMapName;
	182	return *this;
	183	}
	184
	185	//____________________________________________________________________
	186	const char*
	187	AliCentralMultiplicityTask::Manager::GetFullFileName(UShort_t what,
	188	UShort_t sys,
	189	UShort_t sNN,
	190	Short_t field) const
	191	{
	192	return Form("%s/%s",
	193	what == 0 ? GetSecMapPath() : GetAcceptancePath(),
	194	GetFileName(what, sys, sNN, field));
	195	}
	196
	197	//____________________________________________________________________
	198	const char*
	199	AliCentralMultiplicityTask::Manager::GetFileName(UShort_t what ,
	200	UShort_t sys,
	201	UShort_t sNN,
	202	Short_t field) const
	203	{
	204	// Must be static - otherwise the data may disappear on return from
	205	// this member function
	206	static TString fname = "";
6f791cc3	207
6f791cc3	208	switch(what) {
3e478dba	209	case 0: fname.Append(fSecMapName.Data()); break;
3e478dba	210	case 1: fname.Append(fAcceptanceName.Data()); break;
6f791cc3	211	default:
3e478dba	212	::Error("GetFileName",
3e478dba	213	"Invalid indentifier %d for central object, must be 0 or 1!", what);
6f791cc3	214	break;
	215	}
	216	fname.Append(Form("_%s_%04dGeV_%c%1dkG.root",
	217	AliForwardUtil::CollisionSystemString(sys),
	218	sNN, (field < 0 ? 'm' : 'p'), TMath::Abs(field)));
	219
	220	return fname.Data();
6f791cc3	221	}
3e478dba	222
6f791cc3	223	//____________________________________________________________________
3e478dba	224	TH2D*
	225	AliCentralMultiplicityTask::Manager::GetSecMapCorrection(UShort_t vtxbin) const
	226	{
	227	if (!fSecmap) {
	228	::Warning("GetSecMapCorrection","No secondary map defined");
	229	return 0;
	230	}
	231	return fSecmap->GetCorrection(vtxbin);
	232	}
	233	//____________________________________________________________________
	234	TH1D*
	235	AliCentralMultiplicityTask::Manager::GetAcceptanceCorrection(UShort_t vtxbin)
	236	const
	237	{
	238	if (!fAcceptance) {
	239	::Warning("GetAcceptanceCorrection","No acceptance map defined");
	240	return 0;
	241	}
	242	return fAcceptance->GetCorrection(vtxbin);
	243	}
	244
	245	//____________________________________________________________________
	246	void
	247	AliCentralMultiplicityTask::Manager::Init(UShort_t sys,
	248	UShort_t sNN,
	249	Short_t field)
	250	{
6f791cc3	251	TFile fsec(GetFullFileName(0,sys,sNN,field));
3e478dba	252	fSecmap =
3e478dba	253	dynamic_cast<AliCentralCorrSecondaryMap*>(fsec.Get(fSecMapName.Data()));
6f791cc3	254	if(!fSecmap) {
3e478dba	255	::Error("Init", "no central Secondary Map found!") ;
6f791cc3	256	return;
	257	}
	258	TFile facc(GetFullFileName(1,sys,sNN,field));
3e478dba	259	fAcceptance =
	260	dynamic_cast<AliCentralCorrAcceptance*>(facc.Get(fAcceptanceName.Data()));
	261	if(!fAcceptance) {
	262	::Error("Init", "no central Acceptance found!") ;
	263	return;
	264	}
6f791cc3	265
	266	}
	267	//
	268	// EOF
	269	//