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


#include "AliFMDHistCollector.h"
#include <AliESDFMD.h>
#include <TAxis.h>
#include <TList.h>
#include <TMath.h>
#include "AliFMDAnaParameters.h"
#include "AliLog.h"
#include <TH2D.h>
#include <TH1I.h>
#include <TProfile.h>
#include <TProfile2D.h>
#include <TArrayI.h>

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


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

  fNCutBins       = o.fNCutBins;
  fCorrectionCut  = o.fCorrectionCut;
  fFirstBins      = o.fFirstBins;
  fLastBins       = o.fLastBins;
  
  return *this;
}

//____________________________________________________________________
void
AliFMDHistCollector::Init(const TAxis& vtxAxis)
{
  AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();

  Int_t nVz = vtxAxis.GetNbins();
  fFirstBins.Set(5*nVz);
  fLastBins.Set(5*nVz);

  // Find the eta bin ranges 
  for (Int_t iVz = 0; iVz < nVz; iVz++) {
    // Find the first and last eta bin to use for each ring for 
    // each vertex bin.   This is instead of using the methods 
    // provided by AliFMDAnaParameters 
    for (Int_t iIdx = 0; iIdx < 5; iIdx++) {
      UShort_t d = 0;
      Char_t   r = 0;
      GetDetRing(iIdx, d, r);
      
      // Get the background object 
      TH2F* bg    = pars->GetBackgroundCorrection(d,r,iVz);
      Int_t nEta  = bg->GetNbinsX();
      Int_t first = nEta+1;
      Int_t last  = 0;

      // Loop over the eta bins 
      for (Int_t ie = 1; ie <= nEta; ie++) { 
	if (bg->GetBinContent(ie,1) < fCorrectionCut) continue;

	// Loop over the phi bins to make sure that we 
	// do not have holes in the coverage 
	bool ok = true;
	for (Int_t ip = 1; ip <= bg->GetNbinsY(); ip++) { 
	  if (bg->GetBinContent(ie,ip) < 0.001) {
	    ok = false;
	    continue;
	  }
	}
	if (!ok) continue;

	first = TMath::Min(ie, first);
	last  = TMath::Max(ie, last);
      }
      
      // Store the result for later use 
      fFirstBins[iVz*5+iIdx] = first;
      fLastBins[iVz*5+iIdx]  = last;
    } // for j 
  }
}

//____________________________________________________________________
Int_t
AliFMDHistCollector::GetIdx(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;
  }
  return idx;
}
//____________________________________________________________________
void
AliFMDHistCollector::GetDetRing(Int_t idx, UShort_t& d, Char_t& r) const
{
  d = 0; 
  r = '\0';
  switch (idx) { 
  case 0: d = 1; r = 'I'; break;
  case 1: d = 2; r = 'I'; break;
  case 2: d = 2; r = 'O'; break;
  case 3: d = 3; r = 'I'; break;
  case 4: d = 3; r = 'O'; break;
  }
}

//____________________________________________________________________
void
AliFMDHistCollector::GetFirstAndLast(Int_t idx, Int_t vtxbin, 
				     Int_t& first, Int_t& last) const
{
  first = 0; 
  last  = 0;
  
  if (idx < 0) return;
  idx += vtxbin * 5;
      
  if (idx < 0 || idx >= fFirstBins.GetSize()) return;
  
  first = fFirstBins.At(idx)+fNCutBins;  
  last  = fLastBins.At(idx)-fNCutBins;
}

//____________________________________________________________________
Int_t
AliFMDHistCollector::GetFirst(Int_t idx, Int_t v) const 
{
  Int_t f, l;
  GetFirstAndLast(idx,v,f,l);
  return f;
}


//____________________________________________________________________
Int_t
AliFMDHistCollector::GetLast(Int_t idx, Int_t v) const 
{
  Int_t f, l;
  GetFirstAndLast(idx,v,f,l);
  return l;
}

//____________________________________________________________________
Int_t 
AliFMDHistCollector::GetOverlap(UShort_t d, Char_t r, 
				Int_t bin, Int_t vtxbin) const
{
  // Get the possibly overlapping histogram 
  Int_t other = -1;
  if (d == 1) {
    if (bin <= GetLast(2,'I',vtxbin)) other = GetIdx(2,'I');
  }
  else if (d == 2 && r == 'I') {
    if      (bin <= GetLast(2,  'O', vtxbin)) other = GetIdx(2, 'O');
    else if (bin >= GetFirst(1, 'I', vtxbin)) other = GetIdx(1, 'I');
  }
  else if (d == 2 && r == 'O') {
    if (bin >= GetFirst(2, 'I', vtxbin)) other = GetIdx(2,'I');
  }
  else if (d == 3 && r == 'O') {
    if (bin <= GetLast(3, 'I', vtxbin)) other = GetIdx(3, 'I');
  }
  else if (d == 3 && r == 'I') {
    if (bin >= GetFirst(3, 'O', vtxbin)) other = GetIdx(3, 'O');
  }
  // AliInfo(Form("FMD%d%c (%d) -> %d", d, r, GetIdx(d,r), other));
  return other;
}
//____________________________________________________________________
Int_t
AliFMDHistCollector::GetOverlap(Int_t idx, Int_t bin, Int_t vtxbin) const
{
  UShort_t d = 0; 
  Char_t   r = '\0';
  GetDetRing(idx, d, r);
  if (d == 0 || r == '\0') return 0;

  return GetOverlap(d, r, bin, vtxbin);
}
  
  
//____________________________________________________________________
Bool_t
AliFMDHistCollector::Collect(AliForwardUtil::Histos& hists,
			     Int_t                   vtxbin, 
			     TH2D&                   out)
{
  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);
      TH2D*       t = static_cast<TH2D*>(h->Clone(Form("FMD%d%c_tmp",d,r)));
      
      // Outer rings have better phi segmentation - rebin to same as inner. 
      if (q == 1) t->RebinY(2);

      Int_t first = 0;
      Int_t last  = 0;
      GetFirstAndLast(d, r, vtxbin, first, last);

      // Now update profile output 
      for (Int_t iEta = first; iEta <= last; iEta++) { 

	// Get the possibly overlapping histogram 
	Int_t overlap = GetOverlap(d,r,iEta,vtxbin);

	// Fill eta acceptance for this event into the phi underlow bin
	Float_t ooc      = out.GetBinContent(iEta,0);
	Float_t noc      = overlap >= 0? 0.5 : 1;
	out.SetBinContent(iEta, 0, ooc + noc);

	for (Int_t iPhi = 1; iPhi <= h->GetNbinsY(); iPhi++) { 
	  Double_t c = t->GetBinContent(iEta,iPhi);
	  Double_t e = t->GetBinError(iEta,iPhi);

	  // If there's no signal, continue 
	  // if (e <= 0) continue;
	  if (c <= 0 || e <= 0)     continue;
	  
	  // If there's no overlapping histogram (ring), then 
	  // fill in data and continue to the next phi bin 
	  if (overlap < 0) { 
	    out.SetBinContent(iEta,iPhi,c);
	    out.SetBinError(iEta,iPhi,e);
	    continue;
	  }

	  // Get the current bin content and error 
	  Double_t oc = out.GetBinContent(iEta,iPhi);
	  Double_t oe = out.GetBinError(iEta,iPhi);

#define USE_STRAIGHT_MEAN
// #define USE_STRAIGHT_MEAN_NONZERO
// #define USE_WEIGHTED_MEAN
// #define USE_MOST_CERTAIN
#if defined(USE_STRAIGHT_MEAN)
	  // calculate the average of old value (half the original), 
	  // and this value, as well as the summed squared errors 
	  // of the existing content (sqrt((e_1/2)^2=sqrt(e_1^2/4)=e_1/2) 
	  // and half the error of this.   
	  //
	  // So, on the first overlapping histogram we get 
	  // 
	  //    c = c_1 / 2
	  //    e = sqrt((e_1 / 2)^2) = e_1/2
	  // 
	  // On the second we get 
	  // 
	  //    c' = c_2 / 2 + c = c_2 / 2 + c_1 / 2 = (c_1+c_2)/2 
	  //    e' = sqrt(e^2 + (e_2/2)^2) 
	  //       = sqrt(e_1^2/4 + e_2^2/4) 
	  //       = sqrt(1/4 * (e_1^2+e_2^2)) 
	  //       = 1/2 * sqrt(e_1^2 + e_2^2)
	  out.SetBinContent(iEta,iPhi,oc + c/2);
	  out.SetBinError(iEta,iPhi,TMath::Sqrt(oe*oe+(e*e)/4));
#elif defined(USE_STRAIGHT_MEAN_NONZERO)
# define ZERO_OTHER
	  // If there's data in the overlapping histogram, 
	  // calculate the average and add the errors in 
	  // quadrature.  
	  // 
	  // If there's no data in the overlapping histogram, 
	  // then just fill in the data 
	  if (oe > 0) {
	    out.SetBinContent(iEta,iPhi,(oc + c)/2);
	    out.SetBinError(iEta,iPhi,TMath::Sqrt(oe*oe + e*e)/2);
	  }
	  else {
	    out.SetBinContent(iEta,iPhi,c);
	    out.SetBinError(iEta,iPhi,e);
	  }	    
#elif defined(USE_WEIGHTED_MEAN) 
	  // Calculate the weighted mean 
	  Double_t w  = 1/(e*e);
	  Double_t sc = w * c;
	  Double_t sw = w;
	  if (oe > 0) {
	    Double_t ow = 1/(oe*oe);
	    sc          += ow * oc;
	    sw          += ow;
	  }
	  Double_t nc = sc / sw;
	  Double_t ne = TMath::Sqrt(1 / sw);
	  out.SetBinContent(iEta,iPhi,nc,ne);
#elif defined(USE_MOST_CERTAIN)
# define ZERO_OTHER
	  if (e < oe) {
	    out.SetBinContent(iEta,iPhi,c);
	    out.SetBinError(iEta,iPhi,e);
	  }
	  else {
	    out.SetBinContent(iEta,iPhi,oc);
	    out.SetBinError(iEta,iPhi,oe);
	  }
#else 
#         error No method for defining content of overlapping bins defined
#endif
#if defined(ZERO_OTHER)
	  // Get the content of the overlapping histogram, 
	  // and zero the content so that we won't use it 
	  // again 
	  UShort_t od; Char_t oR; 
	  GetDetRing(overlap, od, oR);
	  TH2D* other = hists.Get(od,oR);
	  other->SetBinContent(iEta,iPhi,0);
	  other->SetBinError(iEta,iPhi,0);
#endif
	}
      }
      // Remove temporary histogram 
      delete t;
    } // for r
  } // for d 
  return kTRUE;
}


//____________________________________________________________________
//
// EOF
//
Commit	Line	Data
7e4038b5	1
	2	#include "AliFMDHistCollector.h"
	3	#include <AliESDFMD.h>
	4	#include <TAxis.h>
	5	#include <TList.h>
	6	#include <TMath.h>
	7	#include "AliFMDAnaParameters.h"
	8	#include "AliLog.h"
	9	#include <TH2D.h>
	10	#include <TH1I.h>
	11	#include <TProfile.h>
	12	#include <TProfile2D.h>
	13	#include <TArrayI.h>
	14
	15	ClassImp(AliFMDHistCollector)
	16	#if 0
	17	; // For Emacs
	18	#endif
	19
	20
	21	//____________________________________________________________________
	22	AliFMDHistCollector&
	23	AliFMDHistCollector::operator=(const AliFMDHistCollector& o)
	24	{
	25	SetName(o.GetName());
	26	SetTitle(o.GetTitle());
	27
	28	fNCutBins = o.fNCutBins;
	29	fCorrectionCut = o.fCorrectionCut;
	30	fFirstBins = o.fFirstBins;
	31	fLastBins = o.fLastBins;
	32
	33	return *this;
	34	}
	35
	36	//____________________________________________________________________
	37	void
	38	AliFMDHistCollector::Init(const TAxis& vtxAxis)
	39	{
	40	AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
	41
	42	Int_t nVz = vtxAxis.GetNbins();
	43	fFirstBins.Set(5*nVz);
	44	fLastBins.Set(5*nVz);
	45
	46	// Find the eta bin ranges
	47	for (Int_t iVz = 0; iVz < nVz; iVz++) {
	48	// Find the first and last eta bin to use for each ring for
	49	// each vertex bin. This is instead of using the methods
	50	// provided by AliFMDAnaParameters
	51	for (Int_t iIdx = 0; iIdx < 5; iIdx++) {
	52	UShort_t d = 0;
	53	Char_t r = 0;
	54	GetDetRing(iIdx, d, r);
	55
	56	// Get the background object
	57	TH2F* bg = pars->GetBackgroundCorrection(d,r,iVz);
	58	Int_t nEta = bg->GetNbinsX();
	59	Int_t first = nEta+1;
	60	Int_t last = 0;
	61
	62	// Loop over the eta bins
	63	for (Int_t ie = 1; ie <= nEta; ie++) {
	64	if (bg->GetBinContent(ie,1) < fCorrectionCut) continue;
65
66	// Loop over the phi bins to make sure that we
67	// do not have holes in the coverage
68	bool ok = true;
69	for (Int_t ip = 1; ip <= bg->GetNbinsY(); ip++) {
70	if (bg->GetBinContent(ie,ip) < 0.001) {
71	ok = false;
72	continue;
73	}
74	}
75	if (!ok) continue;
76
77	first = TMath::Min(ie, first);
78	last = TMath::Max(ie, last);
79	}
80
81	// Store the result for later use
82	fFirstBins[iVz*5+iIdx] = first;
83	fLastBins[iVz*5+iIdx] = last;
84	} // for j
85	}
86	}
87
88	//____________________________________________________________________
89	Int_t
90	AliFMDHistCollector::GetIdx(UShort_t d, Char_t r) const
91	{
92	Int_t idx = -1;
93	switch (d) {
94	case 1: idx = 0; break;
95	case 2: idx = 1 + (r == 'I' \|\| r == 'i' ? 0 : 1); break;
96	case 3: idx = 3 + (r == 'I' \|\| r == 'i' ? 0 : 1); break;
97	}
98	return idx;
99	}
100	//____________________________________________________________________
101	void
102	AliFMDHistCollector::GetDetRing(Int_t idx, UShort_t& d, Char_t& r) const
103	{
104	d = 0;
105	r = '\0';
106	switch (idx) {
107	case 0: d = 1; r = 'I'; break;
108	case 1: d = 2; r = 'I'; break;
109	case 2: d = 2; r = 'O'; break;
110	case 3: d = 3; r = 'I'; break;
111	case 4: d = 3; r = 'O'; break;
112	}
113	}
114
115	//____________________________________________________________________
116	void
117	AliFMDHistCollector::GetFirstAndLast(Int_t idx, Int_t vtxbin,
118	Int_t& first, Int_t& last) const
119	{
120	first = 0;
121	last = 0;
122
123	if (idx < 0) return;
124	idx += vtxbin * 5;
125
126	if (idx < 0 \|\| idx >= fFirstBins.GetSize()) return;
127
128	first = fFirstBins.At(idx)+fNCutBins;
129	last = fLastBins.At(idx)-fNCutBins;
130	}
131
132	//____________________________________________________________________
133	Int_t
134	AliFMDHistCollector::GetFirst(Int_t idx, Int_t v) const
135	{
136	Int_t f, l;
137	GetFirstAndLast(idx,v,f,l);
138	return f;
139	}
140
141
142	//____________________________________________________________________
143	Int_t
144	AliFMDHistCollector::GetLast(Int_t idx, Int_t v) const
145	{
146	Int_t f, l;
147	GetFirstAndLast(idx,v,f,l);
148	return l;
149	}
150
151	//____________________________________________________________________
152	Int_t
153	AliFMDHistCollector::GetOverlap(UShort_t d, Char_t r,
154	Int_t bin, Int_t vtxbin) const
155	{
156	// Get the possibly overlapping histogram
157	Int_t other = -1;
158	if (d == 1) {
159	if (bin <= GetLast(2,'I',vtxbin)) other = GetIdx(2,'I');
160	}
161	else if (d == 2 && r == 'I') {
162	if (bin <= GetLast(2, 'O', vtxbin)) other = GetIdx(2, 'O');
163	else if (bin >= GetFirst(1, 'I', vtxbin)) other = GetIdx(1, 'I');
164	}
165	else if (d == 2 && r == 'O') {
166	if (bin >= GetFirst(2, 'I', vtxbin)) other = GetIdx(2,'I');
167	}
168	else if (d == 3 && r == 'O') {
169	if (bin <= GetLast(3, 'I', vtxbin)) other = GetIdx(3, 'I');
170	}
171	else if (d == 3 && r == 'I') {
172	if (bin >= GetFirst(3, 'O', vtxbin)) other = GetIdx(3, 'O');
173	}
174	// AliInfo(Form("FMD%d%c (%d) -> %d", d, r, GetIdx(d,r), other));
175	return other;
176	}
177	//____________________________________________________________________
178	Int_t
179	AliFMDHistCollector::GetOverlap(Int_t idx, Int_t bin, Int_t vtxbin) const
180	{
181	UShort_t d = 0;
182	Char_t r = '\0';
183	GetDetRing(idx, d, r);
184	if (d == 0 \|\| r == '\0') return 0;
185
186	return GetOverlap(d, r, bin, vtxbin);
187	}
188
189
190	//____________________________________________________________________
191	Bool_t
192	AliFMDHistCollector::Collect(AliForwardUtil::Histos& hists,
193	Int_t vtxbin,
194	TH2D& out)
195	{
196	for (UShort_t d=1; d<=3; d++) {
197	UShort_t nr = (d == 1 ? 1 : 2);
198	for (UShort_t q=0; q<nr; q++) {
199	Char_t r = (q == 0 ? 'I' : 'O');
200	TH2D* h = hists.Get(d,r);
201	TH2D* t = static_cast<TH2D*>(h->Clone(Form("FMD%d%c_tmp",d,r)));
202
203	// Outer rings have better phi segmentation - rebin to same as inner.
204	if (q == 1) t->RebinY(2);
205
206	Int_t first = 0;
207	Int_t last = 0;
208	GetFirstAndLast(d, r, vtxbin, first, last);
209
210	// Now update profile output
211	for (Int_t iEta = first; iEta <= last; iEta++) {
212
213	// Get the possibly overlapping histogram
214	Int_t overlap = GetOverlap(d,r,iEta,vtxbin);
215
216	// Fill eta acceptance for this event into the phi underlow bin
217	Float_t ooc = out.GetBinContent(iEta,0);
218	Float_t noc = overlap >= 0? 0.5 : 1;
219	out.SetBinContent(iEta, 0, ooc + noc);
220
221	for (Int_t iPhi = 1; iPhi <= h->GetNbinsY(); iPhi++) {
222	Double_t c = t->GetBinContent(iEta,iPhi);
223	Double_t e = t->GetBinError(iEta,iPhi);
224
225	// If there's no signal, continue
226	// if (e <= 0) continue;
227	if (c <= 0 \|\| e <= 0) continue;
228
229	// If there's no overlapping histogram (ring), then
230	// fill in data and continue to the next phi bin
231	if (overlap < 0) {
232	out.SetBinContent(iEta,iPhi,c);
233	out.SetBinError(iEta,iPhi,e);
234	continue;
235	}
236
237	// Get the current bin content and error
238	Double_t oc = out.GetBinContent(iEta,iPhi);
239	Double_t oe = out.GetBinError(iEta,iPhi);
240
241	#define USE_STRAIGHT_MEAN
242	// #define USE_STRAIGHT_MEAN_NONZERO
243	// #define USE_WEIGHTED_MEAN
244	// #define USE_MOST_CERTAIN
245	#if defined(USE_STRAIGHT_MEAN)
246	// calculate the average of old value (half the original),
247	// and this value, as well as the summed squared errors
248	// of the existing content (sqrt((e_1/2)^2=sqrt(e_1^2/4)=e_1/2)
249	// and half the error of this.
250	//
251	// So, on the first overlapping histogram we get
252	//
253	// c = c_1 / 2
254	// e = sqrt((e_1 / 2)^2) = e_1/2
255	//
256	// On the second we get
257	//
258	// c' = c_2 / 2 + c = c_2 / 2 + c_1 / 2 = (c_1+c_2)/2
259	// e' = sqrt(e^2 + (e_2/2)^2)
260	// = sqrt(e_1^2/4 + e_2^2/4)
261	// = sqrt(1/4 * (e_1^2+e_2^2))
262	// = 1/2 * sqrt(e_1^2 + e_2^2)
263	out.SetBinContent(iEta,iPhi,oc + c/2);
264	out.SetBinError(iEta,iPhi,TMath::Sqrt(oeoe+(ee)/4));
265	#elif defined(USE_STRAIGHT_MEAN_NONZERO)
266	# define ZERO_OTHER
267	// If there's data in the overlapping histogram,
268	// calculate the average and add the errors in
269	// quadrature.
270	//
271	// If there's no data in the overlapping histogram,
272	// then just fill in the data
273	if (oe > 0) {
274	out.SetBinContent(iEta,iPhi,(oc + c)/2);
275	out.SetBinError(iEta,iPhi,TMath::Sqrt(oeoe + ee)/2);
276	}
277	else {
278	out.SetBinContent(iEta,iPhi,c);
279	out.SetBinError(iEta,iPhi,e);
280	}
281	#elif defined(USE_WEIGHTED_MEAN)
282	// Calculate the weighted mean
283	Double_t w = 1/(e*e);
284	Double_t sc = w * c;
285	Double_t sw = w;
286	if (oe > 0) {
287	Double_t ow = 1/(oe*oe);
288	sc += ow * oc;
289	sw += ow;
290	}
291	Double_t nc = sc / sw;
292	Double_t ne = TMath::Sqrt(1 / sw);
293	out.SetBinContent(iEta,iPhi,nc,ne);
294	#elif defined(USE_MOST_CERTAIN)
295	# define ZERO_OTHER
296	if (e < oe) {
297	out.SetBinContent(iEta,iPhi,c);
298	out.SetBinError(iEta,iPhi,e);
299	}
300	else {
301	out.SetBinContent(iEta,iPhi,oc);
302	out.SetBinError(iEta,iPhi,oe);
303	}
304	#else
305	# error No method for defining content of overlapping bins defined
306	#endif
307	#if defined(ZERO_OTHER)
308	// Get the content of the overlapping histogram,
309	// and zero the content so that we won't use it
310	// again
311	UShort_t od; Char_t oR;
312	GetDetRing(overlap, od, oR);
313	TH2D* other = hists.Get(od,oR);
314	other->SetBinContent(iEta,iPhi,0);
315	other->SetBinError(iEta,iPhi,0);
316	#endif
317	}
318	}
319	// Remove temporary histogram
320	delete t;
321	} // for r
322	} // for d
323	return kTRUE;
324	}
325
326
327	//____________________________________________________________________
328	//
329	// EOF
330	//
331
332
333