[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 <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;
  Int_t vtxbin = (Int_t)vertexBinDouble +1 ; //HHD: The vtxbins are 1-10, not 0-9
  
  // 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 error = aodnew*TMath::Sqrt(TMath::Power(aodHist->GetBinError(nx,ny)/aodvalue,2) +
					 TMath::Power(hAcceptance->GetBinError(nx)/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")
{


}
//____________________________________________________________________
const char* AliCentralMultiplicityTask::Manager::GetFileName(UShort_t  what ,
							     UShort_t  sys, 
							     UShort_t  sNN,
							     Short_t   field) {
  
  TString fname = "";
  
  switch(what) {
  case 0:
    fname.Append(fSecMapName.Data());
    break;
  case 1: 
    fname.Append(fAcceptanceName.Data());
    break;
  default:
    printf("Invalid indentifier for central object, must be 0 or 1!");
    break;
  }
  fname.Append(Form("_%s_%04dGeV_%c%1dkG.root", 
		    AliForwardUtil::CollisionSystemString(sys), 
		    sNN, (field < 0 ? 'm' : 'p'), TMath::Abs(field)));
  
  return fname.Data();
  
}
//____________________________________________________________________
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) {
    printf("no central Secondary Map found!") ;
    return;
  }
  TFile facc(GetFullFileName(1,sys,sNN,field));
  fAcceptance = dynamic_cast<AliCentralCorrAcceptance*>(facc.Get(fAcceptanceName.Data()));
if(!fAcceptance) {
  printf("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>
	27	#include <iostream>
	28	#include <iomanip>
	29
	30	//====================================================================
	31	AliCentralMultiplicityTask::AliCentralMultiplicityTask(const char* name)
	32	: AliAnalysisTaskSE(name),
	33	fData(0),
	34	fList(0),
	35	fAODCentral(kFALSE),
	36	fManager()
	37	{
	38
	39	DefineOutput(1, TList::Class());
	40	}
	41	//____________________________________________________________________
	42	void AliCentralMultiplicityTask::UserCreateOutputObjects() {
	43
	44	AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
	45	AliAODHandler* ah =
	46	dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
	47	if (!ah) AliFatal("No AOD output handler set in analysis manager");
	48
	49
	50	TObject* obj = &fAODCentral;
	51	ah->AddBranch("AliAODCentralMult", &obj);
	52
	53	fList = new TList();
	54	PostData(1,fList);
	55
	56	}
	57	//____________________________________________________________________
	58	void AliCentralMultiplicityTask::UserExec(Option_t* /option/) {
	59
	60	AliESDInputHandler* eventHandler =
	61	dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()
	62	->GetInputEventHandler());
	63	if (!eventHandler) {
	64	AliWarning("No inputhandler found for this event!");
65	return;}
66
67	AliESDEvent* esd = eventHandler->GetEvent();
68
69	//Selecting only events with \|valid vertex\| < 10 cm
70	const AliESDVertex* vertex = esd->GetPrimaryVertexSPD();
71	if (!vertex) return;
72	if(!(vertex->GetStatus())) return;
73	if(vertex->GetNContributors() <= 0) return ;
74	if(vertex->GetZRes() > 0.1 ) return;
75	Double_t vertexXYZ[3]={0,0,0};
76	vertex->GetXYZ(vertexXYZ);
77	if(TMath::Abs(vertexXYZ[2]) > 10) return;
78
79	Double_t delta = 2 ;
80	Double_t vertexBinDouble = (vertexXYZ[2] + 10) / delta;
81	Int_t vtxbin = (Int_t)vertexBinDouble +1 ; //HHD: The vtxbins are 1-10, not 0-9
82
83	// Make sure AOD is filled
84	AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
85	AliAODHandler* ah =
86	dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
87	if (!ah)
88	AliFatal("No AOD output handler set in analysis manager");
89
90	ah->SetFillAOD(kTRUE);
91
92	//Doing analysis
93	fAODCentral.Clear("");
94	TH2D *aodHist = &(fAODCentral.GetHistogram());
95
96	const AliMultiplicity* spdmult = esd->GetMultiplicity();
97	//Filling clusters in layer 1 used for tracklets...
98	for(Int_t j = 0; j< spdmult->GetNumberOfTracklets();j++) {
99	aodHist->Fill(spdmult->GetEta(j),spdmult->GetPhi(j));
100	}
101	//...and then the unused ones in layer 1
102	for(Int_t j = 0; j< spdmult->GetNumberOfSingleClusters();j++) {
103	aodHist->Fill(-TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.)),
104	spdmult->GetPhiSingle(j));
105	}
106
107	// Corrections
108
109	TH2D* hSecMap = fManager.GetSecMapCorrection(vtxbin);
110	TH1D* hAcceptance = fManager.GetAcceptanceCorrection(vtxbin);
111
112	aodHist->Divide(hSecMap);
113
114	for(Int_t nx = 1; nx <= aodHist->GetNbinsX(); nx++) {
115	Float_t acccor = hAcceptance->GetBinContent(nx);
116
117	Bool_t etabinSeen = kFALSE;
118	for(Int_t ny = 1; ny <= aodHist->GetNbinsY(); ny++) {
119	Float_t aodvalue = aodHist->GetBinContent(nx,ny);
120	Float_t seccor = hSecMap->GetBinContent(nx,ny);
121	if(seccor > 0.5) etabinSeen = kTRUE;
122	if(aodvalue < 0.000001) { aodHist->SetBinContent(nx,ny, 0); continue; }
123
124	Float_t aodnew = aodvalue / acccor ;
125	aodHist->SetBinContent(nx,ny, aodnew);
126	Float_t error = aodnew*TMath::Sqrt(TMath::Power(aodHist->GetBinError(nx,ny)/aodvalue,2) +
127	TMath::Power(hAcceptance->GetBinError(nx)/acccor,2) );
128	aodHist->SetBinError(nx,ny,error);
129
130	}
131	//Filling underflow bin if we eta bin is in range
132	if(etabinSeen) aodHist->SetBinContent(nx,0, 1.);
133
134	}
135
136	PostData(1,fList);
137
138	}
139	//____________________________________________________________________
140	void AliCentralMultiplicityTask::Terminate(Option_t* /option/) {
141
142
143
144	}
145	//____________________________________________________________________
146	void
147	AliCentralMultiplicityTask::Print(Option_t* /option/) const
148	{
149
150	}
151	//____________________________________________________________________
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	//____________________________________________________________________
164	const char* AliCentralMultiplicityTask::Manager::GetFileName(UShort_t what ,
165	UShort_t sys,
166	UShort_t sNN,
167	Short_t field) {
168
169	TString fname = "";
170
171	switch(what) {
172	case 0:
173	fname.Append(fSecMapName.Data());
174	break;
175	case 1:
176	fname.Append(fAcceptanceName.Data());
177	break;
178	default:
179	printf("Invalid indentifier for central object, must be 0 or 1!");
180	break;
181	}
182	fname.Append(Form("_%s_%04dGeV_%c%1dkG.root",
183	AliForwardUtil::CollisionSystemString(sys),
184	sNN, (field < 0 ? 'm' : 'p'), TMath::Abs(field)));
185
186	return fname.Data();
187
188	}
189	//____________________________________________________________________
190	void AliCentralMultiplicityTask::Manager::Init( UShort_t sys,
191	UShort_t sNN,
192	Short_t field) {
193
194	TFile fsec(GetFullFileName(0,sys,sNN,field));
195	fSecmap = dynamic_cast<AliCentralCorrSecondaryMap*>(fsec.Get(fSecMapName.Data()));
196	if(!fSecmap) {
197	printf("no central Secondary Map found!") ;
198	return;
199	}
200	TFile facc(GetFullFileName(1,sys,sNN,field));
201	fAcceptance = dynamic_cast<AliCentralCorrAcceptance*>(facc.Get(fAcceptanceName.Data()));
202	if(!fAcceptance) {
203	printf("no central Acceptance found!") ;
204	return;
205	}
206
207	}
208	//
209	// EOF
210	//