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

#include "AliFMDCorrections.h"
#include <AliESDFMD.h>
#include <TAxis.h>
#include <TList.h>
#include <TMath.h>
#include "AliFMDAnaParameters.h"
#include "AliLog.h"
#include <TH2D.h>

ClassImp(AliFMDCorrections)
#if 0
; // For Emacs
#endif 

//____________________________________________________________________
AliFMDCorrections::AliFMDCorrections()
  : TNamed(), 
    fRingHistos(),
    fMultCut(0.3)
{}

//____________________________________________________________________
AliFMDCorrections::AliFMDCorrections(const char* title)
  : TNamed("fmdCorrections", title), 
    fRingHistos(), 
    fMultCut(0.3)
{
  fRingHistos.SetName(GetName());
  fRingHistos.Add(new RingHistos(1, 'I'));
  fRingHistos.Add(new RingHistos(2, 'I'));
  fRingHistos.Add(new RingHistos(2, 'O'));
  fRingHistos.Add(new RingHistos(3, 'I'));
  fRingHistos.Add(new RingHistos(3, 'O'));
}

//____________________________________________________________________
AliFMDCorrections::AliFMDCorrections(const AliFMDCorrections& o)
  : TNamed(o), 
    fRingHistos(), 
    fMultCut(o.fMultCut)
{
  TIter    next(&o.fRingHistos);
  TObject* obj = 0;
  while ((obj = next())) fRingHistos.Add(obj);
}

//____________________________________________________________________
AliFMDCorrections::~AliFMDCorrections()
{
  fRingHistos.Delete();
}

//____________________________________________________________________
AliFMDCorrections&
AliFMDCorrections::operator=(const AliFMDCorrections& o)
{
  SetName(o.GetName());
  SetTitle(o.GetTitle());

  fMultCut = o.fMultCut;

  fRingHistos.Delete();
  TIter    next(&o.fRingHistos);
  TObject* obj = 0;
  while ((obj = next())) fRingHistos.Add(obj);
  
  return *this;
}

//____________________________________________________________________
AliFMDCorrections::RingHistos*
AliFMDCorrections::GetRingHistos(UShort_t d, Char_t r) const
{
  Int_t idx = -1;
  switch (d) { 
  case 1: idx = 0; break;
  case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
  case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
  }
  if (idx < 0 || idx >= fRingHistos.GetEntries()) return 0;
  
  return static_cast<RingHistos*>(fRingHistos.At(idx));
}
    
//____________________________________________________________________
Bool_t
AliFMDCorrections::Correct(AliForwardUtil::Histos& hists,
			   Int_t                   vtxbin)
{
  AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();

  for (UShort_t d=1; d<=3; d++) { 
    UShort_t nr = (d == 1 ? 1 : 2);
    for (UShort_t q=0; q<nr; q++) { 
      Char_t      r = (q == 0 ? 'I' : 'O');
      TH2D*       h = hists.Get(d,r);
      RingHistos* rh= GetRingHistos(d,r);
      TH2F*       bg= pars->GetBackgroundCorrection(d, r, vtxbin);
      TH2F*       ef= pars->GetEventSelectionEfficiency("INEL",vtxbin,r);
      if (!bg) { 
	AliWarning(Form("No secondary correction for FMDM%d%c in vertex bin %d",
			d, r, vtxbin));
	continue;
      }
      if (!ef) { 
	AliWarning(Form("No event %s selection efficiency in vertex bin %d",
			"INEL", vtxbin));
	continue;
      }

      // Divide by primary/total ratio
      h->Divide(bg);
      
      // Divide by the event selection efficiency 
      h->Divide(ef);

      if(!pars->SharingEffPresent()) { 
	AliWarning("No sharing efficiencies");
	continue;
      }

      TH1F* sf = pars->GetSharingEfficiencyTrVtx(d,r,vtxbin); 
      
      for (Int_t ieta = 1; ieta <= h->GetNbinsX(); ieta++) {
	Float_t c  = sf->GetBinContent(ieta);
	Float_t ec = sf->GetBinError(ieta);
	
	if (c == 0) continue;

	for (Int_t iphi = 1; iphi <= h->GetNbinsY(); iphi++) { 
	  Double_t m  = h->GetBinContent(ieta, iphi) / c;
	  Double_t em = h->GetBinError(ieta, iphi);
	  
	  Double_t e  = TMath::Sqrt(em * em + (m * ec) * (m * ec));
	  
	  h->SetBinContent(ieta,iphi,m);
	  h->SetBinError(ieta,iphi,e);
	}
      }

      rh->fDensity->Add(h);
    }
  }
  
  return kTRUE;
}

//____________________________________________________________________
void
AliFMDCorrections::ScaleHistograms(TList* dir, Int_t nEvents)
{
  if (nEvents <= 0) return;
  TList* d = static_cast<TList*>(dir->FindObject(GetName()));
  if (!d) return;

  TIter    next(&fRingHistos);
  RingHistos* o = 0;
  while ((o = static_cast<RingHistos*>(next())))
    o->ScaleHistograms(d, nEvents);
}
//____________________________________________________________________
void
AliFMDCorrections::DefineOutput(TList* dir)
{
  TList* d = new TList;
  d->SetName(GetName());
  dir->Add(d);
  TIter    next(&fRingHistos);
  RingHistos* o = 0;
  while ((o = static_cast<RingHistos*>(next()))) {
    o->Output(d);
  }
}

//====================================================================
AliFMDCorrections::RingHistos::RingHistos()
  : AliForwardUtil::RingHistos(), 
    fDensity(0)
{}

//____________________________________________________________________
AliFMDCorrections::RingHistos::RingHistos(UShort_t d, Char_t r)
  : AliForwardUtil::RingHistos(d,r), 
    fDensity(0)
{
  fDensity = new TH2D(Form("FMD%d%c_Primary_Density", d, r), 
		      Form("in FMD%d%c", d, r), 
		      200, -4, 6, (r == 'I' || r == 'i' ? 20 : 40), 
		      0, 2*TMath::Pi());
  fDensity->SetDirectory(0);
  fDensity->SetXTitle("#eta");
  fDensity->SetYTitle("#phi [radians]");
  fDensity->SetZTitle("Primary N_{ch} density");
}
//____________________________________________________________________
AliFMDCorrections::RingHistos::RingHistos(const RingHistos& o)
  : AliForwardUtil::RingHistos(o), 
    fDensity(o.fDensity)
{}

//____________________________________________________________________
AliFMDCorrections::RingHistos&
AliFMDCorrections::RingHistos::operator=(const RingHistos& o)
{
  AliForwardUtil::RingHistos::operator=(o);
  
  if (fDensity) delete fDensity;
  
  fDensity = static_cast<TH2D*>(o.fDensity->Clone());
  
  return *this;
}
//____________________________________________________________________
AliFMDCorrections::RingHistos::~RingHistos()
{
  if (fDensity) delete fDensity;
}

//____________________________________________________________________
void
AliFMDCorrections::RingHistos::Output(TList* dir)
{
  TList* d = DefineOutputList(dir);
  d->Add(fDensity);
}

//____________________________________________________________________
void
AliFMDCorrections::RingHistos::ScaleHistograms(TList* dir, Int_t nEvents)
{ 
  TList* l = GetOutputList(dir);
  if (!l) return; 

  TH1* density = GetOutputHist(l,Form("%s_Primary_Density", fName.Data()));
  if (density) density->Scale(1./nEvents);
}

//____________________________________________________________________
//
// EOF
//
Commit	Line	Data
7e4038b5	1	#include "AliFMDCorrections.h"
	2	#include <AliESDFMD.h>
	3	#include <TAxis.h>
	4	#include <TList.h>
	5	#include <TMath.h>
	6	#include "AliFMDAnaParameters.h"
	7	#include "AliLog.h"
	8	#include <TH2D.h>
	9
	10	ClassImp(AliFMDCorrections)
	11	#if 0
	12	; // For Emacs
	13	#endif
	14
	15	//____________________________________________________________________
	16	AliFMDCorrections::AliFMDCorrections()
	17	: TNamed(),
	18	fRingHistos(),
	19	fMultCut(0.3)
	20	{}
	21
	22	//____________________________________________________________________
	23	AliFMDCorrections::AliFMDCorrections(const char* title)
	24	: TNamed("fmdCorrections", title),
	25	fRingHistos(),
	26	fMultCut(0.3)
	27	{
	28	fRingHistos.SetName(GetName());
	29	fRingHistos.Add(new RingHistos(1, 'I'));
	30	fRingHistos.Add(new RingHistos(2, 'I'));
	31	fRingHistos.Add(new RingHistos(2, 'O'));
	32	fRingHistos.Add(new RingHistos(3, 'I'));
	33	fRingHistos.Add(new RingHistos(3, 'O'));
	34	}
	35
	36	//____________________________________________________________________
	37	AliFMDCorrections::AliFMDCorrections(const AliFMDCorrections& o)
	38	: TNamed(o),
	39	fRingHistos(),
	40	fMultCut(o.fMultCut)
	41	{
	42	TIter next(&o.fRingHistos);
	43	TObject* obj = 0;
	44	while ((obj = next())) fRingHistos.Add(obj);
	45	}
	46
	47	//____________________________________________________________________
	48	AliFMDCorrections::~AliFMDCorrections()
	49	{
	50	fRingHistos.Delete();
	51	}
	52
	53	//____________________________________________________________________
	54	AliFMDCorrections&
	55	AliFMDCorrections::operator=(const AliFMDCorrections& o)
	56	{
	57	SetName(o.GetName());
	58	SetTitle(o.GetTitle());
	59
	60	fMultCut = o.fMultCut;
	61
	62	fRingHistos.Delete();
	63	TIter next(&o.fRingHistos);
	64	TObject* obj = 0;
65	while ((obj = next())) fRingHistos.Add(obj);
66
67	return *this;
68	}
69
70	//____________________________________________________________________
71	AliFMDCorrections::RingHistos*
72	AliFMDCorrections::GetRingHistos(UShort_t d, Char_t r) const
73	{
74	Int_t idx = -1;
75	switch (d) {
76	case 1: idx = 0; break;
77	case 2: idx = 1 + (r == 'I' \|\| r == 'i' ? 0 : 1); break;
78	case 3: idx = 3 + (r == 'I' \|\| r == 'i' ? 0 : 1); break;
79	}
80	if (idx < 0 \|\| idx >= fRingHistos.GetEntries()) return 0;
81
82	return static_cast<RingHistos*>(fRingHistos.At(idx));
83	}
84
85	//____________________________________________________________________
86	Bool_t
87	AliFMDCorrections::Correct(AliForwardUtil::Histos& hists,
88	Int_t vtxbin)
89	{
90	AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
91
92	for (UShort_t d=1; d<=3; d++) {
93	UShort_t nr = (d == 1 ? 1 : 2);
94	for (UShort_t q=0; q<nr; q++) {
95	Char_t r = (q == 0 ? 'I' : 'O');
96	TH2D* h = hists.Get(d,r);
97	RingHistos* rh= GetRingHistos(d,r);
98	TH2F* bg= pars->GetBackgroundCorrection(d, r, vtxbin);
99	TH2F* ef= pars->GetEventSelectionEfficiency("INEL",vtxbin,r);
100	if (!bg) {
101	AliWarning(Form("No secondary correction for FMDM%d%c in vertex bin %d",
102	d, r, vtxbin));
103	continue;
104	}
105	if (!ef) {
106	AliWarning(Form("No event %s selection efficiency in vertex bin %d",
107	"INEL", vtxbin));
108	continue;
109	}
110
111	// Divide by primary/total ratio
112	h->Divide(bg);
113
114	// Divide by the event selection efficiency
115	h->Divide(ef);
116
117	if(!pars->SharingEffPresent()) {
118	AliWarning("No sharing efficiencies");
119	continue;
120	}
121
122	TH1F* sf = pars->GetSharingEfficiencyTrVtx(d,r,vtxbin);
123
124	for (Int_t ieta = 1; ieta <= h->GetNbinsX(); ieta++) {
125	Float_t c = sf->GetBinContent(ieta);
126	Float_t ec = sf->GetBinError(ieta);
127
128	if (c == 0) continue;
129
130	for (Int_t iphi = 1; iphi <= h->GetNbinsY(); iphi++) {
131	Double_t m = h->GetBinContent(ieta, iphi) / c;
132	Double_t em = h->GetBinError(ieta, iphi);
133
134	Double_t e = TMath::Sqrt(em * em + (m * ec) * (m * ec));
135
136	h->SetBinContent(ieta,iphi,m);
137	h->SetBinError(ieta,iphi,e);
138	}
139	}
140
141	rh->fDensity->Add(h);
142	}
143	}
144
145	return kTRUE;
146	}
147
148	//____________________________________________________________________
149	void
9d99b0dd	150	AliFMDCorrections::ScaleHistograms(TList* dir, Int_t nEvents)
7e4038b5	151	{
7e4038b5	152	if (nEvents <= 0) return;
9d99b0dd	153	TList* d = static_cast<TList*>(dir->FindObject(GetName()));
9d99b0dd	154	if (!d) return;
7e4038b5	155
	156	TIter next(&fRingHistos);
	157	RingHistos* o = 0;
9d99b0dd	158	while ((o = static_cast<RingHistos*>(next())))
9d99b0dd	159	o->ScaleHistograms(d, nEvents);
7e4038b5	160	}
7e4038b5	161	//____________________________________________________________________
7e4038b5	162	void
9d99b0dd	163	AliFMDCorrections::DefineOutput(TList* dir)
7e4038b5	164	{
	165	TList* d = new TList;
	166	d->SetName(GetName());
	167	dir->Add(d);
	168	TIter next(&fRingHistos);
	169	RingHistos* o = 0;
	170	while ((o = static_cast<RingHistos*>(next()))) {
	171	o->Output(d);
	172	}
	173	}
	174
	175	//====================================================================
	176	AliFMDCorrections::RingHistos::RingHistos()
9d99b0dd	177	: AliForwardUtil::RingHistos(),
7e4038b5	178	fDensity(0)
	179	{}
	180
	181	//____________________________________________________________________
	182	AliFMDCorrections::RingHistos::RingHistos(UShort_t d, Char_t r)
9d99b0dd	183	: AliForwardUtil::RingHistos(d,r),
7e4038b5	184	fDensity(0)
	185	{
	186	fDensity = new TH2D(Form("FMD%d%c_Primary_Density", d, r),
	187	Form("in FMD%d%c", d, r),
	188	200, -4, 6, (r == 'I' \|\| r == 'i' ? 20 : 40),
	189	0, 2*TMath::Pi());
	190	fDensity->SetDirectory(0);
	191	fDensity->SetXTitle("#eta");
	192	fDensity->SetYTitle("#phi [radians]");
	193	fDensity->SetZTitle("Primary N_{ch} density");
	194	}
	195	//____________________________________________________________________
	196	AliFMDCorrections::RingHistos::RingHistos(const RingHistos& o)
9d99b0dd	197	: AliForwardUtil::RingHistos(o),
7e4038b5	198	fDensity(o.fDensity)
	199	{}
	200
	201	//____________________________________________________________________
	202	AliFMDCorrections::RingHistos&
	203	AliFMDCorrections::RingHistos::operator=(const RingHistos& o)
	204	{
9d99b0dd	205	AliForwardUtil::RingHistos::operator=(o);
7e4038b5	206
	207	if (fDensity) delete fDensity;
	208
	209	fDensity = static_cast<TH2D*>(o.fDensity->Clone());
	210
	211	return *this;
	212	}
	213	//____________________________________________________________________
	214	AliFMDCorrections::RingHistos::~RingHistos()
	215	{
	216	if (fDensity) delete fDensity;
	217	}
	218
	219	//____________________________________________________________________
	220	void
	221	AliFMDCorrections::RingHistos::Output(TList* dir)
	222	{
9d99b0dd	223	TList* d = DefineOutputList(dir);
7e4038b5	224	d->Add(fDensity);
9d99b0dd	225	}
	226
	227	//____________________________________________________________________
	228	void
	229	AliFMDCorrections::RingHistos::ScaleHistograms(TList* dir, Int_t nEvents)
	230	{
	231	TList* l = GetOutputList(dir);
	232	if (!l) return;
	233
	234	TH1* density = GetOutputHist(l,Form("%s_Primary_Density", fName.Data()));
	235	if (density) density->Scale(1./nEvents);
7e4038b5	236	}
	237
	238	//____________________________________________________________________
	239	//
	240	// EOF
	241	//
	242
	243
	244