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

// 
// This class collects the event histograms into single histograms, 
// one for each ring in each vertex bin.  
//
// Input:
//   - AliESDFMD object possibly corrected for sharing
//
// Output:
//   - 5 RingHistos objects - each with a number of vertex dependent 
//     2D histograms of the inclusive charge particle density 
// 
// HistCollector used: 
//   - AliFMDCorrSecondaryMap
//
#include "AliFMDHistCollector.h"
#include <AliESDFMD.h>
#include <TAxis.h>
#include <TList.h>
#include <TMath.h>
#include "AliForwardCorrectionManager.h"
#include "AliLog.h"
#include <TH2D.h>
#include <TH1I.h>
#include <TProfile.h>
#include <TProfile2D.h>
#include <TArrayI.h>
#include <TROOT.h>
#include <iostream>
#include <iomanip>

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


//____________________________________________________________________
AliFMDHistCollector&
AliFMDHistCollector::operator=(const AliFMDHistCollector& o)
{
  // 
  // Assignement operator
  // 
  // Parameters:
  //    o Object to assign from 
  //
  // Return:
  //    Reference to this object
  //
  TNamed::operator=(o);

  fNCutBins       = o.fNCutBins;
  fCorrectionCut  = o.fCorrectionCut;
  fFirstBins      = o.fFirstBins;
  fLastBins       = o.fLastBins;
  fDebug          = o.fDebug;

  return *this;
}

//____________________________________________________________________
void
AliFMDHistCollector::Init(const TAxis& vtxAxis)
{
  // 
  // Intialise 
  // 
  // Parameters:
  //    vtxAxis  Vertex axis 
  //  

  AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();

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

  // Find the eta bin ranges 
  for (UShort_t iVz = 1; 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);
      TH2D* bg    = fcm.GetSecondaryMap()->GetCorrection(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-1)*5+iIdx] = first;
      fLastBins[(iVz-1)*5+iIdx]  = last;
    } // for j 
  }
}

//____________________________________________________________________
Int_t
AliFMDHistCollector::GetIdx(UShort_t d, Char_t r) const
{
  // 
  // Get the ring index from detector number and ring identifier 
  // 
  // Parameters:
  //    d Detector
  //    r Ring identifier 
  // 
  // Return:
  //    ring index or -1 in case of problems 
  //
  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
{
  // 
  // Get the detector and ring from the ring index 
  // 
  // Parameters:
  //    idx Ring index 
  //    d   On return, the detector or 0 in case of errors 
  //    r   On return, the ring id or '0' in case of errors 
  //
  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, UShort_t vtxbin, 
				     Int_t& first, Int_t& last) const
{
  // 
  // Get the first and last eta bin to use for a given ring and vertex 
  // 
  // Parameters:
  //    idx      Ring index as given by GetIdx
  //    vtxBin   Vertex bin (1 based) 
  //    first    On return, the first eta bin to use 
  //    last     On return, the last eta bin to use 
  //
  first = 0; 
  last  = 0;
  
  if (idx    <  0) return;
  if (vtxbin <= 0) return;
  idx += (vtxbin-1) * 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, UShort_t v) const 
{
  // 
  // Get the first eta bin to use for a given ring and vertex 
  // 
  // Parameters:
  //    idx Ring index as given by GetIdx
  //    v vertex bin (1 based)
  // 
  // Return:
  //    First eta bin to use, or -1 in case of problems 
  //  
  Int_t f, l;
  GetFirstAndLast(idx,v,f,l);
  return f;
}


//____________________________________________________________________
Int_t
AliFMDHistCollector::GetLast(Int_t idx, UShort_t v) const 
{
  // 
  // Get the last eta bin to use for a given ring and vertex 
  // 
  // Parameters:
  //    idx Ring index as given by GetIdx
  //    v vertex bin (1 based)
  // 
  // Return:
  //    Last eta bin to use, or -1 in case of problems 
  //  
  Int_t f, l;
  GetFirstAndLast(idx,v,f,l);
  return l;
}

//____________________________________________________________________
Int_t 
AliFMDHistCollector::GetOverlap(UShort_t d, Char_t r, 
				Int_t bin,  UShort_t vtxbin) const
{
  // 
  // Get the possibly overlapping histogram of eta bin @a e in 
  // detector and ring 
  // 
  // Parameters:
  //    d Detector
  //    r Ring 
  //    e Eta bin
  //    v Vertex bin (1 based)
  //
  // Return:
  //    Overlapping histogram index or -1
  //

  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, UShort_t vtxbin) const
{
  // 
  // Get the possibly overlapping histogram of eta bin @a e in 
  // detector and ring 
  // 
  // Parameters:
  //    i Ring index
  //    e Eta bin
  //    v Vertex bin (1 based)
  //
  // Return:
  //    Overlapping histogram index or -1
  //
  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,
			     UShort_t                vtxbin, 
			     TH2D&                   out)
{
  // 
  // Do the calculations 
  // 
  // Parameters:
  //    hists    Cache of histograms 
  //    vtxBin   Vertex bin (1 based)
  //    out      Output histogram
  // 
  // Return:
  //    true on successs 
  //
  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;
}

//____________________________________________________________________
void
AliFMDHistCollector::Print(Option_t* /* option */) const
{
  // 
  // Print information 
  // 
  // Parameters:
  //    option Not used
  //
  char ind[gROOT->GetDirLevel()+1];
  for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
  ind[gROOT->GetDirLevel()] = '\0';
  std::cout << ind << "AliFMDHistCollector: " << GetName() << '\n'
	    << ind << " # of cut bins:          " << fNCutBins << '\n'
	    << ind << " Correction cut:         " << fCorrectionCut << '\n'
	    << ind << " Bin ranges:\n" << ind << "  v_z bin";
  Int_t nVz = fFirstBins.fN / 5;
  for (UShort_t iVz = 1; iVz <= nVz; iVz++) 
    std::cout << " | " << std::setw(7) << iVz;
  std::cout << '\n' << ind << "  / ring ";
  for (UShort_t iVz = 1; iVz <= nVz; iVz++) 
    std::cout << "-+--------";
  std::cout << std::endl;
    
  for (Int_t iIdx = 0; iIdx < 5; iIdx++) {
    UShort_t d = 0;
    Char_t   r = 0;
    GetDetRing(iIdx, d, r);
    
    std::cout << ind << "  FMD" << d << r << "  ";
    for (UShort_t iVz = 1; iVz <= nVz; iVz++) {
      Int_t first, last;
      GetFirstAndLast(iIdx, iVz, first, last);
      std::cout << " | " << std::setw(3) << first << "-" 
		<< std::setw(3) << last;
    }
    std::cout << std::endl;
  }
}

//____________________________________________________________________
//
// EOF
//
Commit	Line	Data
7984e5f7	1	//
	2	// This class collects the event histograms into single histograms,
	3	// one for each ring in each vertex bin.
	4	//
	5	// Input:
	6	// - AliESDFMD object possibly corrected for sharing
	7	//
	8	// Output:
	9	// - 5 RingHistos objects - each with a number of vertex dependent
	10	// 2D histograms of the inclusive charge particle density
	11	//
	12	// HistCollector used:
	13	// - AliFMDCorrSecondaryMap
	14	//
7e4038b5	15	#include "AliFMDHistCollector.h"
	16	#include <AliESDFMD.h>
	17	#include <TAxis.h>
	18	#include <TList.h>
	19	#include <TMath.h>
0bd4b00f	20	#include "AliForwardCorrectionManager.h"
7e4038b5	21	#include "AliLog.h"
	22	#include <TH2D.h>
	23	#include <TH1I.h>
	24	#include <TProfile.h>
	25	#include <TProfile2D.h>
	26	#include <TArrayI.h>
0bd4b00f	27	#include <TROOT.h>
	28	#include <iostream>
	29	#include <iomanip>
7e4038b5	30
	31	ClassImp(AliFMDHistCollector)
	32	#if 0
	33	; // For Emacs
	34	#endif
	35
	36
	37	//____________________________________________________________________
	38	AliFMDHistCollector&
	39	AliFMDHistCollector::operator=(const AliFMDHistCollector& o)
	40	{
7984e5f7	41	//
	42	// Assignement operator
	43	//
	44	// Parameters:
	45	// o Object to assign from
	46	//
	47	// Return:
	48	// Reference to this object
	49	//
ea3e5d95	50	TNamed::operator=(o);
7e4038b5	51
	52	fNCutBins = o.fNCutBins;
	53	fCorrectionCut = o.fCorrectionCut;
	54	fFirstBins = o.fFirstBins;
	55	fLastBins = o.fLastBins;
ea3e5d95	56	fDebug = o.fDebug;
ea3e5d95	57
7e4038b5	58	return *this;
	59	}
	60
	61	//____________________________________________________________________
	62	void
	63	AliFMDHistCollector::Init(const TAxis& vtxAxis)
	64	{
7984e5f7	65	//
	66	// Intialise
	67	//
	68	// Parameters:
	69	// vtxAxis Vertex axis
	70	//
	71
0bd4b00f	72	AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
7e4038b5	73
0bd4b00f	74	UShort_t nVz = vtxAxis.GetNbins();
7e4038b5	75	fFirstBins.Set(5*nVz);
	76	fLastBins.Set(5*nVz);
	77
	78	// Find the eta bin ranges
0bd4b00f	79	for (UShort_t iVz = 1; iVz <= nVz; iVz++) {
7e4038b5	80	// Find the first and last eta bin to use for each ring for
	81	// each vertex bin. This is instead of using the methods
	82	// provided by AliFMDAnaParameters
	83	for (Int_t iIdx = 0; iIdx < 5; iIdx++) {
	84	UShort_t d = 0;
	85	Char_t r = 0;
	86	GetDetRing(iIdx, d, r);
	87
	88	// Get the background object
0bd4b00f	89	// TH2F* bg = pars->GetBackgroundCorrection(d,r,iVz);
0bd4b00f	90	TH2D* bg = fcm.GetSecondaryMap()->GetCorrection(d,r,iVz);
7e4038b5	91	Int_t nEta = bg->GetNbinsX();
	92	Int_t first = nEta+1;
	93	Int_t last = 0;
	94
	95	// Loop over the eta bins
	96	for (Int_t ie = 1; ie <= nEta; ie++) {
	97	if (bg->GetBinContent(ie,1) < fCorrectionCut) continue;
	98
	99	// Loop over the phi bins to make sure that we
	100	// do not have holes in the coverage
	101	bool ok = true;
	102	for (Int_t ip = 1; ip <= bg->GetNbinsY(); ip++) {
	103	if (bg->GetBinContent(ie,ip) < 0.001) {
	104	ok = false;
	105	continue;
	106	}
	107	}
	108	if (!ok) continue;
	109
	110	first = TMath::Min(ie, first);
	111	last = TMath::Max(ie, last);
	112	}
	113
	114	// Store the result for later use
0bd4b00f	115	fFirstBins[(iVz-1)*5+iIdx] = first;
0bd4b00f	116	fLastBins[(iVz-1)*5+iIdx] = last;
7e4038b5	117	} // for j
	118	}
	119	}
	120
	121	//____________________________________________________________________
	122	Int_t
	123	AliFMDHistCollector::GetIdx(UShort_t d, Char_t r) const
	124	{
7984e5f7	125	//
	126	// Get the ring index from detector number and ring identifier
	127	//
	128	// Parameters:
	129	// d Detector
	130	// r Ring identifier
	131	//
	132	// Return:
	133	// ring index or -1 in case of problems
	134	//
7e4038b5	135	Int_t idx = -1;
	136	switch (d) {
	137	case 1: idx = 0; break;
	138	case 2: idx = 1 + (r == 'I' \|\| r == 'i' ? 0 : 1); break;
	139	case 3: idx = 3 + (r == 'I' \|\| r == 'i' ? 0 : 1); break;
	140	}
	141	return idx;
	142	}
	143	//____________________________________________________________________
	144	void
	145	AliFMDHistCollector::GetDetRing(Int_t idx, UShort_t& d, Char_t& r) const
	146	{
7984e5f7	147	//
	148	// Get the detector and ring from the ring index
	149	//
	150	// Parameters:
	151	// idx Ring index
	152	// d On return, the detector or 0 in case of errors
	153	// r On return, the ring id or '0' in case of errors
	154	//
7e4038b5	155	d = 0;
	156	r = '\0';
	157	switch (idx) {
	158	case 0: d = 1; r = 'I'; break;
	159	case 1: d = 2; r = 'I'; break;
	160	case 2: d = 2; r = 'O'; break;
	161	case 3: d = 3; r = 'I'; break;
	162	case 4: d = 3; r = 'O'; break;
	163	}
	164	}
	165
	166	//____________________________________________________________________
	167	void
0bd4b00f	168	AliFMDHistCollector::GetFirstAndLast(Int_t idx, UShort_t vtxbin,
7e4038b5	169	Int_t& first, Int_t& last) const
7e4038b5	170	{
7984e5f7	171	//
	172	// Get the first and last eta bin to use for a given ring and vertex
	173	//
	174	// Parameters:
	175	// idx Ring index as given by GetIdx
	176	// vtxBin Vertex bin (1 based)
	177	// first On return, the first eta bin to use
	178	// last On return, the last eta bin to use
	179	//
7e4038b5	180	first = 0;
	181	last = 0;
	182
0bd4b00f	183	if (idx < 0) return;
	184	if (vtxbin <= 0) return;
	185	idx += (vtxbin-1) * 5;
7e4038b5	186
	187	if (idx < 0 \|\| idx >= fFirstBins.GetSize()) return;
	188
	189	first = fFirstBins.At(idx)+fNCutBins;
	190	last = fLastBins.At(idx)-fNCutBins;
	191	}
	192
	193	//____________________________________________________________________
	194	Int_t
0bd4b00f	195	AliFMDHistCollector::GetFirst(Int_t idx, UShort_t v) const
7e4038b5	196	{
7984e5f7	197	//
	198	// Get the first eta bin to use for a given ring and vertex
	199	//
	200	// Parameters:
	201	// idx Ring index as given by GetIdx
	202	// v vertex bin (1 based)
	203	//
	204	// Return:
	205	// First eta bin to use, or -1 in case of problems
	206	//
7e4038b5	207	Int_t f, l;
	208	GetFirstAndLast(idx,v,f,l);
	209	return f;
	210	}
	211
	212
	213	//____________________________________________________________________
	214	Int_t
0bd4b00f	215	AliFMDHistCollector::GetLast(Int_t idx, UShort_t v) const
7e4038b5	216	{
7984e5f7	217	//
	218	// Get the last eta bin to use for a given ring and vertex
	219	//
	220	// Parameters:
	221	// idx Ring index as given by GetIdx
	222	// v vertex bin (1 based)
	223	//
	224	// Return:
	225	// Last eta bin to use, or -1 in case of problems
	226	//
7e4038b5	227	Int_t f, l;
	228	GetFirstAndLast(idx,v,f,l);
	229	return l;
	230	}
	231
	232	//____________________________________________________________________
	233	Int_t
	234	AliFMDHistCollector::GetOverlap(UShort_t d, Char_t r,
0bd4b00f	235	Int_t bin, UShort_t vtxbin) const
7e4038b5	236	{
7984e5f7	237	//
	238	// Get the possibly overlapping histogram of eta bin @a e in
	239	// detector and ring
	240	//
	241	// Parameters:
	242	// d Detector
	243	// r Ring
	244	// e Eta bin
	245	// v Vertex bin (1 based)
	246	//
	247	// Return:
	248	// Overlapping histogram index or -1
	249	//
	250
7e4038b5	251	Int_t other = -1;
	252	if (d == 1) {
	253	if (bin <= GetLast(2,'I',vtxbin)) other = GetIdx(2,'I');
	254	}
	255	else if (d == 2 && r == 'I') {
	256	if (bin <= GetLast(2, 'O', vtxbin)) other = GetIdx(2, 'O');
	257	else if (bin >= GetFirst(1, 'I', vtxbin)) other = GetIdx(1, 'I');
	258	}
	259	else if (d == 2 && r == 'O') {
	260	if (bin >= GetFirst(2, 'I', vtxbin)) other = GetIdx(2,'I');
	261	}
	262	else if (d == 3 && r == 'O') {
	263	if (bin <= GetLast(3, 'I', vtxbin)) other = GetIdx(3, 'I');
	264	}
	265	else if (d == 3 && r == 'I') {
	266	if (bin >= GetFirst(3, 'O', vtxbin)) other = GetIdx(3, 'O');
	267	}
	268	// AliInfo(Form("FMD%d%c (%d) -> %d", d, r, GetIdx(d,r), other));
	269	return other;
	270	}
	271	//____________________________________________________________________
	272	Int_t
0bd4b00f	273	AliFMDHistCollector::GetOverlap(Int_t idx, Int_t bin, UShort_t vtxbin) const
7e4038b5	274	{
7984e5f7	275	//
	276	// Get the possibly overlapping histogram of eta bin @a e in
	277	// detector and ring
	278	//
	279	// Parameters:
	280	// i Ring index
	281	// e Eta bin
	282	// v Vertex bin (1 based)
	283	//
	284	// Return:
	285	// Overlapping histogram index or -1
	286	//
7e4038b5	287	UShort_t d = 0;
	288	Char_t r = '\0';
	289	GetDetRing(idx, d, r);
	290	if (d == 0 \|\| r == '\0') return 0;
	291
	292	return GetOverlap(d, r, bin, vtxbin);
	293	}
	294
	295
	296	//____________________________________________________________________
	297	Bool_t
	298	AliFMDHistCollector::Collect(AliForwardUtil::Histos& hists,
0bd4b00f	299	UShort_t vtxbin,
7e4038b5	300	TH2D& out)
7e4038b5	301	{
7984e5f7	302	//
	303	// Do the calculations
	304	//
	305	// Parameters:
	306	// hists Cache of histograms
	307	// vtxBin Vertex bin (1 based)
	308	// out Output histogram
	309	//
	310	// Return:
	311	// true on successs
	312	//
7e4038b5	313	for (UShort_t d=1; d<=3; d++) {
	314	UShort_t nr = (d == 1 ? 1 : 2);
	315	for (UShort_t q=0; q<nr; q++) {
	316	Char_t r = (q == 0 ? 'I' : 'O');
	317	TH2D* h = hists.Get(d,r);
	318	TH2D* t = static_cast<TH2D*>(h->Clone(Form("FMD%d%c_tmp",d,r)));
	319
	320	// Outer rings have better phi segmentation - rebin to same as inner.
	321	if (q == 1) t->RebinY(2);
	322
	323	Int_t first = 0;
	324	Int_t last = 0;
	325	GetFirstAndLast(d, r, vtxbin, first, last);
	326
	327	// Now update profile output
	328	for (Int_t iEta = first; iEta <= last; iEta++) {
	329
	330	// Get the possibly overlapping histogram
	331	Int_t overlap = GetOverlap(d,r,iEta,vtxbin);
	332
	333	// Fill eta acceptance for this event into the phi underlow bin
	334	Float_t ooc = out.GetBinContent(iEta,0);
	335	Float_t noc = overlap >= 0? 0.5 : 1;
	336	out.SetBinContent(iEta, 0, ooc + noc);
	337
	338	for (Int_t iPhi = 1; iPhi <= h->GetNbinsY(); iPhi++) {
	339	Double_t c = t->GetBinContent(iEta,iPhi);
	340	Double_t e = t->GetBinError(iEta,iPhi);
	341
	342	// If there's no signal, continue
	343	// if (e <= 0) continue;
	344	if (c <= 0 \|\| e <= 0) continue;
	345
	346	// If there's no overlapping histogram (ring), then
	347	// fill in data and continue to the next phi bin
	348	if (overlap < 0) {
	349	out.SetBinContent(iEta,iPhi,c);
	350	out.SetBinError(iEta,iPhi,e);
	351	continue;
	352	}
	353
	354	// Get the current bin content and error
	355	Double_t oc = out.GetBinContent(iEta,iPhi);
	356	Double_t oe = out.GetBinError(iEta,iPhi);
	357
	358	#define USE_STRAIGHT_MEAN
	359	// #define USE_STRAIGHT_MEAN_NONZERO
	360	// #define USE_WEIGHTED_MEAN
	361	// #define USE_MOST_CERTAIN
	362	#if defined(USE_STRAIGHT_MEAN)
	363	// calculate the average of old value (half the original),
	364	// and this value, as well as the summed squared errors
	365	// of the existing content (sqrt((e_1/2)^2=sqrt(e_1^2/4)=e_1/2)
	366	// and half the error of this.
	367	//
	368	// So, on the first overlapping histogram we get
	369	//
	370	// c = c_1 / 2
	371	// e = sqrt((e_1 / 2)^2) = e_1/2
	372	//
	373	// On the second we get
	374	//
	375	// c' = c_2 / 2 + c = c_2 / 2 + c_1 / 2 = (c_1+c_2)/2
	376	// e' = sqrt(e^2 + (e_2/2)^2)
377	// = sqrt(e_1^2/4 + e_2^2/4)
378	// = sqrt(1/4 * (e_1^2+e_2^2))
379	// = 1/2 * sqrt(e_1^2 + e_2^2)
380	out.SetBinContent(iEta,iPhi,oc + c/2);
381	out.SetBinError(iEta,iPhi,TMath::Sqrt(oeoe+(ee)/4));
382	#elif defined(USE_STRAIGHT_MEAN_NONZERO)
383	# define ZERO_OTHER
384	// If there's data in the overlapping histogram,
385	// calculate the average and add the errors in
386	// quadrature.
387	//
388	// If there's no data in the overlapping histogram,
389	// then just fill in the data
390	if (oe > 0) {
391	out.SetBinContent(iEta,iPhi,(oc + c)/2);
392	out.SetBinError(iEta,iPhi,TMath::Sqrt(oeoe + ee)/2);
393	}
394	else {
395	out.SetBinContent(iEta,iPhi,c);
396	out.SetBinError(iEta,iPhi,e);
397	}
398	#elif defined(USE_WEIGHTED_MEAN)
399	// Calculate the weighted mean
400	Double_t w = 1/(e*e);
401	Double_t sc = w * c;
402	Double_t sw = w;
403	if (oe > 0) {
404	Double_t ow = 1/(oe*oe);
405	sc += ow * oc;
406	sw += ow;
407	}
408	Double_t nc = sc / sw;
409	Double_t ne = TMath::Sqrt(1 / sw);
410	out.SetBinContent(iEta,iPhi,nc,ne);
411	#elif defined(USE_MOST_CERTAIN)
412	# define ZERO_OTHER
413	if (e < oe) {
414	out.SetBinContent(iEta,iPhi,c);
415	out.SetBinError(iEta,iPhi,e);
416	}
417	else {
418	out.SetBinContent(iEta,iPhi,oc);
419	out.SetBinError(iEta,iPhi,oe);
420	}
421	#else
422	# error No method for defining content of overlapping bins defined
423	#endif
424	#if defined(ZERO_OTHER)
425	// Get the content of the overlapping histogram,
426	// and zero the content so that we won't use it
427	// again
428	UShort_t od; Char_t oR;
429	GetDetRing(overlap, od, oR);
430	TH2D* other = hists.Get(od,oR);
431	other->SetBinContent(iEta,iPhi,0);
432	other->SetBinError(iEta,iPhi,0);
433	#endif
434	}
435	}
436	// Remove temporary histogram
437	delete t;
438	} // for r
439	} // for d
440	return kTRUE;
441	}
442
0bd4b00f	443	//____________________________________________________________________
	444	void
	445	AliFMDHistCollector::Print(Option_t* /* option */) const
	446	{
7984e5f7	447	//
	448	// Print information
	449	//
	450	// Parameters:
	451	// option Not used
	452	//
0bd4b00f	453	char ind[gROOT->GetDirLevel()+1];
	454	for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
	455	ind[gROOT->GetDirLevel()] = '\0';
	456	std::cout << ind << "AliFMDHistCollector: " << GetName() << '\n'
	457	<< ind << " # of cut bins: " << fNCutBins << '\n'
	458	<< ind << " Correction cut: " << fCorrectionCut << '\n'
	459	<< ind << " Bin ranges:\n" << ind << " v_z bin";
	460	Int_t nVz = fFirstBins.fN / 5;
	461	for (UShort_t iVz = 1; iVz <= nVz; iVz++)
	462	std::cout << " \| " << std::setw(7) << iVz;
	463	std::cout << '\n' << ind << " / ring ";
	464	for (UShort_t iVz = 1; iVz <= nVz; iVz++)
	465	std::cout << "-+--------";
	466	std::cout << std::endl;
	467
	468	for (Int_t iIdx = 0; iIdx < 5; iIdx++) {
	469	UShort_t d = 0;
	470	Char_t r = 0;
	471	GetDetRing(iIdx, d, r);
	472
	473	std::cout << ind << " FMD" << d << r << " ";
	474	for (UShort_t iVz = 1; iVz <= nVz; iVz++) {
	475	Int_t first, last;
	476	GetFirstAndLast(iIdx, iVz, first, last);
	477	std::cout << " \| " << std::setw(3) << first << "-"
	478	<< std::setw(3) << last;
	479	}
	480	std::cout << std::endl;
	481	}
	482	}
7e4038b5	483
	484	//____________________________________________________________________
	485	//
	486	// EOF
	487	//
	488
	489
	490