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

// Calculate the total energy loss 
#include "AliFMDHistCollector.h"
#include <AliESDFMD.h>
#include <TAxis.h>
#include <TList.h>
#include <TMath.h>
// #include "AliFMDAnaParameters.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)
{
  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)
{
  // AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  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
{
  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, UShort_t vtxbin, 
				     Int_t& first, Int_t& last) const
{
  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 
{
  Int_t f, l;
  GetFirstAndLast(idx,v,f,l);
  return f;
}


//____________________________________________________________________
Int_t
AliFMDHistCollector::GetLast(Int_t idx, UShort_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,  UShort_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, UShort_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,
			     UShort_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;
}

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