[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliFMDCorrector.cxx

//
// This class calculates the exclusive charged particle density
// in each for the 5 FMD rings. 
//
#include "AliFMDCorrector.h"
#include <AliESDFMD.h>
#include <TAxis.h>
#include <TList.h>
#include <TMath.h>
#include "AliForwardCorrectionManager.h"
#include "AliFMDCorrSecondaryMap.h"
#include "AliFMDCorrVertexBias.h"
#include "AliFMDCorrMergingEfficiency.h"
#include "AliFMDCorrAcceptance.h"
#include "AliLog.h"
#include <TH2D.h>
#include <TROOT.h>
#include <THStack.h>
#include <iostream>
#include <iomanip>

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

//____________________________________________________________________
AliFMDCorrector::AliFMDCorrector()
  : TNamed(), 
    fRingHistos(),
    fUseSecondaryMap(true),
    fUseVertexBias(false),
    fUseAcceptance(true),
    fUseMergingEfficiency(false),
    fDebug(0)
{
  // Constructor
  DGUARD(fDebug, 3, "Default CTOR of AliFMDCorrector");
}

//____________________________________________________________________
AliFMDCorrector::AliFMDCorrector(const char* title)
  : TNamed("fmdCorrector", title), 
    fRingHistos(), 
    fUseSecondaryMap(true),
    fUseVertexBias(false),
    fUseAcceptance(true),
    fUseMergingEfficiency(false),
    fDebug(0)
{
  // Constructor 
  // 
  // Parameters: 
  //   title   Title
  DGUARD(fDebug, 3, "Named CTOR of AliFMDCorrector: %s", title);
  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'));
}

//____________________________________________________________________
AliFMDCorrector::AliFMDCorrector(const AliFMDCorrector& o)
  : TNamed(o), 
    fRingHistos(), 
    fUseSecondaryMap(o.fUseSecondaryMap),
    fUseVertexBias(o.fUseVertexBias),
    fUseAcceptance(o.fUseAcceptance),
    fUseMergingEfficiency(o.fUseMergingEfficiency),
    fDebug(o.fDebug)
{
  // Copy constructor 
  // 
  // Parameters: 
  //  o  Object to copy from 
  DGUARD(fDebug, 3, "Copy CTOR of AliFMDCorrector");
  TIter    next(&o.fRingHistos);
  TObject* obj = 0;
  while ((obj = next())) fRingHistos.Add(obj);
}

//____________________________________________________________________
AliFMDCorrector::~AliFMDCorrector()
{
  // Destructor 
  // 
  // 
  DGUARD(fDebug, 3, "DTOR of AliFMDCorrector");
  fRingHistos.Delete();
}

//____________________________________________________________________
AliFMDCorrector&
AliFMDCorrector::operator=(const AliFMDCorrector& o)
{
  // Assignment operator 
  // 
  // Parameters:
  //   o   Object to assign from 
  DGUARD(fDebug, 3, "Assignment of AliFMDCorrector");
  if (&o == this) return *this; 
  TNamed::operator=(o);

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

//____________________________________________________________________
void
AliFMDCorrector::SetupForData(const TAxis&)
{
  //
  // Initialize this object
  //
  // Parameters:
  //    etaAxis Eta axis to use
  //
  DGUARD(fDebug, 1, "Initialization of AliFMDCorrector");
  if (!fUseSecondaryMap)
    AliWarning("Secondary maps not used - BE CAREFUL");
  if (!fUseVertexBias)
    AliWarning("Vertex bias not used");
  if (!fUseAcceptance)
    AliWarning("Acceptance from dead-channels not used");
}

//____________________________________________________________________
AliFMDCorrector::RingHistos*
AliFMDCorrector::GetRingHistos(UShort_t d, Char_t r) const
{
  // 
  // Get the ring histogram container 
  // Parameters:
  //    d Detector
  //    r Ring 
  // 
  // Return:
  //    Ring histogram container 
  //
  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));
}
    
//____________________________________________________________________
void
AliFMDCorrector::DivideMap(TH2* num, const TH2* denom,
			   Bool_t alsoUnderOver) const
{
  // 
  // Implement TH1::Divide but 
  // - Assume compatible histograms 
  // - Unless third argument is true, do not divide over/under flow bins
  // 
  if (!num || !denom) return;

  Int_t first = (alsoUnderOver ? 0 : 1);
  Int_t lastX = num->GetNbinsX() + (alsoUnderOver ? 1 : 0);
  Int_t lastY = num->GetNbinsY() + (alsoUnderOver ? 1 : 0);
  
  for (Int_t ix = first; ix <= lastX; ix++) {
    for (Int_t iy = first; iy <= lastY; iy++) { 
      Int_t    bin = num->GetBin(ix,iy);
      Double_t c0  = num->GetBinContent(bin);
      Double_t c1  = denom->GetBinContent(bin);
      if (!c1) { 
	num->SetBinContent(bin,0);
	num->SetBinError(bin, 0);
	continue;
      }
      Double_t w   = c0 / c1;
      Double_t e0  = num->GetBinError(bin);
      Double_t e1  = denom->GetBinError(bin);
      Double_t c12 = c1*c1;
      Double_t e2  = (e0*e0*c1*c1 + e1*e1*c0*c0)/(c12*c12);
      
      num->SetBinContent(bin, w);
      num->SetBinError(bin, TMath::Sqrt(e2));
    }
  }
}
//____________________________________________________________________
Bool_t
AliFMDCorrector::Correct(AliForwardUtil::Histos& hists,
			 UShort_t                vtxbin)
{
  // 
  // Do the calculations 
  // Parameters:
  //    hists    Cache of histograms 
  //    vtxBin   Vertex bin 
  // 
  // Return:
  //    true on successs 
  //
  DGUARD(fDebug, 1, "Correct histograms of AliFMDCorrector");
  AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();

  UShort_t uvb = vtxbin;
  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);

      if (fUseSecondaryMap) {
        TH2D*  bg = fcm.GetSecondaryMap()->GetCorrection(d, r, uvb);
        if (!bg) {
          AliWarning(Form("No secondary correction for FMDM%d%c "
			  "in vertex bin %d", d, r, uvb));
          continue;
        }
        // Divide by primary/total ratio
	DivideMap(h, bg, false);
      }
      if (fUseVertexBias) {
        TH2D*  ef = fcm.GetVertexBias()->GetCorrection(r, uvb);
        if (!ef) {
          AliWarning(Form("No event %s vertex bias correction in vertex bin %d",
                          (r == 'I' || r == 'i' ? "inner" : "outer"), uvb));
          continue;
        }
        // Divide by the event selection efficiency
	DivideMap(h, ef, false);
      }
      if (fUseAcceptance) {
        TH2D*  ac = fcm.GetAcceptance()->GetCorrection(d, r, uvb);
        if (!ac) {
          AliWarning(Form("No acceptance correction for FMD%d%c in "
			  "vertex bin %d", d, r, uvb));
          continue;
        }
	// Fill overflow bin with ones 
	for (Int_t i = 1; i <= h->GetNbinsX(); i++) 
	  h->SetBinContent(i, h->GetNbinsY()+1, 1);

        // Divide by the acceptance correction - 
	DivideMap(h, ac, fcm.GetAcceptance()->HasOverflow());
      }

      if (fUseMergingEfficiency) {
	if (!fcm.GetMergingEfficiency()) { 
	  AliWarning("No merging efficiencies");
	  continue;
	}
	TH1D* sf = fcm.GetMergingEfficiency()->GetCorrection(d,r,uvb);
	if (!fcm.GetMergingEfficiency()->GetCorrection(d,r,uvb)) { 
	  AliWarning(Form("No merging efficiency for FMD%d%c at vertex bin %d",
			  d, r, uvb));
	  continue;
	}
	
      
	for (Int_t ieta = 1; ieta <= h->GetNbinsX(); ieta++) {
	  Float_t c  = sf->GetBinContent(ieta);
	  Float_t ec = sf->GetBinError(ieta);
	  	  
	  for (Int_t iphi = 1; iphi <= h->GetNbinsY(); iphi++) { 
	    if (c == 0) {
	      h->SetBinContent(ieta,iphi,0);
	      h->SetBinError(ieta,iphi,0);
	      continue;
	    }

	    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
AliFMDCorrector::Terminate(const TList* dir, TList* output, Int_t nEvents)
{
  // 
  // Scale the histograms to the total number of events 
  // Parameters:
  //    dir     Where the output is stored
  //    nEvents Number of events 
  //
  DGUARD(fDebug, 1, "Scale histograms of AliFMDCorrector");
  if (nEvents <= 0) return;
  TList* d = static_cast<TList*>(dir->FindObject(GetName()));
  if (!d) return;

  TList* out = new TList;
  out->SetName(d->GetName());
  out->SetOwner();

  TIter    next(&fRingHistos);
  RingHistos* o = 0;
  THStack* sums = new THStack("sums", "Sums of ring results");
  while ((o = static_cast<RingHistos*>(next()))) {
    o->Terminate(d, nEvents);
    if (!o->fDensity) { 
      Warning("Terminate", "No density from %s", o->GetName());
      continue;
    }
    TH1D* sum = o->fDensity->ProjectionX(o->GetName(), 1, 
					 o->fDensity->GetNbinsY(),"e");
    sum->Scale(1., "width");
    sum->SetTitle(o->GetName());
    sum->SetDirectory(0);
    sum->SetYTitle("#sum N_{ch,primary}");
    sums->Add(sum);
  }
  out->Add(sums);
  output->Add(out);
}
//____________________________________________________________________
void
AliFMDCorrector::CreateOutputObjects(TList* dir)
{
  // 
  // Output diagnostic histograms to directory 
  // 
  // Parameters:
  //    dir List to write in
  //  
  DGUARD(fDebug, 1, "Define output of AliFMDCorrector");
  TList* d = new TList;
  d->SetName(GetName());
  dir->Add(d);

  d->Add(AliForwardUtil::MakeParameter("secondary", fUseSecondaryMap));
  d->Add(AliForwardUtil::MakeParameter("vertexBias", fUseVertexBias));
  d->Add(AliForwardUtil::MakeParameter("acceptance", fUseAcceptance));
  d->Add(AliForwardUtil::MakeParameter("merging", fUseMergingEfficiency));
  

  TIter    next(&fRingHistos);
  RingHistos* o = 0;
  while ((o = static_cast<RingHistos*>(next()))) {
    o->CreateOutputObjects(d);
  }
}

#define PF(N,V,...)					\
  AliForwardUtil::PrintField(N,V, ## __VA_ARGS__)
#define PFB(N,FLAG)				\
  do {									\
    AliForwardUtil::PrintName(N);					\
    std::cout << std::boolalpha << (FLAG) << std::noboolalpha << std::endl; \
  } while(false)
#define PFV(N,VALUE)					\
  do {							\
    AliForwardUtil::PrintName(N);			\
    std::cout << (VALUE) << std::endl; } while(false)

//____________________________________________________________________
void
AliFMDCorrector::Print(Option_t* /* option */) const
{
  // 
  // Print information
  // Parameters:
  //    option Not used 
  //
  AliForwardUtil::PrintTask(*this);
  gROOT->IncreaseDirLevel();
  PFB("Use secondary maps",		fUseSecondaryMap );
  PFB("Use vertex bias",		fUseVertexBias );
  PFB("Use acceptance",			fUseAcceptance );
  PFB("Use merging efficiency",		fUseMergingEfficiency);
  gROOT->DecreaseDirLevel();  
}

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

//____________________________________________________________________
AliFMDCorrector::RingHistos::RingHistos(UShort_t d, Char_t r)
  : AliForwardUtil::RingHistos(d,r), 
    fDensity(0)
{
  // 
  // Constructor
  // Parameters:
  //    d detector
  //    r ring 
  //
  fDensity = new TH2D("primaryDensity", 
		      "#sum N_{ch,primary}/(#Delta#eta#Delta#phi)", 
		      200, -4, 6, (r == 'I' || r == 'i' ? 20 : 40), 
		      0, 2*TMath::Pi());
  fDensity->SetDirectory(0);
  fDensity->SetMarkerColor(Color());
  fDensity->Sumw2();
  fDensity->SetXTitle("#eta");
  fDensity->SetYTitle("#phi [radians]");
  fDensity->SetZTitle("Primary N_{ch} density");
}
//____________________________________________________________________
AliFMDCorrector::RingHistos::RingHistos(const RingHistos& o)
  : AliForwardUtil::RingHistos(o), 
    fDensity(o.fDensity)
{
  // 
  // Copy constructor 
  // Parameters:
  //    o Object to copy from 
  //
}

//____________________________________________________________________
AliFMDCorrector::RingHistos&
AliFMDCorrector::RingHistos::operator=(const RingHistos& o)
{
  // 
  // Assignment operator 
  // Parameters:
  //    o Object to assign from 
  // 
  // Return:
  //    Reference to this 
  //
  if (&o == this) return *this; 
  AliForwardUtil::RingHistos::operator=(o);
  
  if (fDensity) delete fDensity;
  
  fDensity = static_cast<TH2D*>(o.fDensity->Clone());
  
  return *this;
}
//____________________________________________________________________
AliFMDCorrector::RingHistos::~RingHistos()
{
  // 
  // Destructor 
  //
  // if (fDensity) delete fDensity;
}

//____________________________________________________________________
void
AliFMDCorrector::RingHistos::CreateOutputObjects(TList* dir)
{
  // 
  // Make output 
  // Parameters:
  //    dir Where to put it 
  //
  TList* d = DefineOutputList(dir);
  d->Add(fDensity);
}

//____________________________________________________________________
void
AliFMDCorrector::RingHistos::Terminate(TList* dir, Int_t nEvents)
{ 
  // 
  // Scale the histograms to the total number of events 
  // Parameters:
  //    dir     where the output is stored
  //    nEvents Number of events 
  //
  TList* l = GetOutputList(dir);
  if (!l) return; 

  TH2D* density = static_cast<TH2D*>(GetOutputHist(l,"primaryDensity"));
  if (density) density->Scale(1./nEvents);
  fDensity = density;
}

//____________________________________________________________________
//
// EOF
//
Commit	Line	Data
7984e5f7	1	//
	2	// This class calculates the exclusive charged particle density
	3	// in each for the 5 FMD rings.
	4	//
72cc12cd	5	#include "AliFMDCorrector.h"
7e4038b5	6	#include <AliESDFMD.h>
	7	#include <TAxis.h>
	8	#include <TList.h>
	9	#include <TMath.h>
0bd4b00f	10	#include "AliForwardCorrectionManager.h"
8449e3e0	11	#include "AliFMDCorrSecondaryMap.h"
fb3430ac	12	#include "AliFMDCorrVertexBias.h"
	13	#include "AliFMDCorrMergingEfficiency.h"
	14	#include "AliFMDCorrAcceptance.h"
7e4038b5	15	#include "AliLog.h"
7e4038b5	16	#include <TH2D.h>
0bd4b00f	17	#include <TROOT.h>
5bb5d1f6	18	#include <THStack.h>
0bd4b00f	19	#include <iostream>
0bd4b00f	20	#include <iomanip>
7e4038b5	21
72cc12cd	22	ClassImp(AliFMDCorrector)
7e4038b5	23	#if 0
	24	; // For Emacs
	25	#endif
	26
	27	//____________________________________________________________________
72cc12cd	28	AliFMDCorrector::AliFMDCorrector()
7e4038b5	29	: TNamed(),
7e4038b5	30	fRingHistos(),
7ec4d843	31	fUseSecondaryMap(true),
9d05ffeb	32	fUseVertexBias(false),
7ec4d843	33	fUseAcceptance(true),
9d05ffeb	34	fUseMergingEfficiency(false),
ea3e5d95	35	fDebug(0)
7984e5f7	36	{
7984e5f7	37	// Constructor
5ca83fee	38	DGUARD(fDebug, 3, "Default CTOR of AliFMDCorrector");
7984e5f7	39	}
7e4038b5	40
7e4038b5	41	//____________________________________________________________________
72cc12cd	42	AliFMDCorrector::AliFMDCorrector(const char* title)
cc83fca2	43	: TNamed("fmdCorrector", title),
7e4038b5	44	fRingHistos(),
7ec4d843	45	fUseSecondaryMap(true),
9d05ffeb	46	fUseVertexBias(false),
7ec4d843	47	fUseAcceptance(true),
9d05ffeb	48	fUseMergingEfficiency(false),
ea3e5d95	49	fDebug(0)
7e4038b5	50	{
7984e5f7	51	// Constructor
	52	//
	53	// Parameters:
	54	// title Title
5ca83fee	55	DGUARD(fDebug, 3, "Named CTOR of AliFMDCorrector: %s", title);
7e4038b5	56	fRingHistos.SetName(GetName());
	57	fRingHistos.Add(new RingHistos(1, 'I'));
	58	fRingHistos.Add(new RingHistos(2, 'I'));
	59	fRingHistos.Add(new RingHistos(2, 'O'));
	60	fRingHistos.Add(new RingHistos(3, 'I'));
	61	fRingHistos.Add(new RingHistos(3, 'O'));
	62	}
	63
	64	//____________________________________________________________________
72cc12cd	65	AliFMDCorrector::AliFMDCorrector(const AliFMDCorrector& o)
7e4038b5	66	: TNamed(o),
7e4038b5	67	fRingHistos(),
7ec4d843	68	fUseSecondaryMap(o.fUseSecondaryMap),
	69	fUseVertexBias(o.fUseVertexBias),
	70	fUseAcceptance(o.fUseAcceptance),
81eda625	71	fUseMergingEfficiency(o.fUseMergingEfficiency),
ea3e5d95	72	fDebug(o.fDebug)
7e4038b5	73	{
7984e5f7	74	// Copy constructor
	75	//
	76	// Parameters:
	77	// o Object to copy from
5ca83fee	78	DGUARD(fDebug, 3, "Copy CTOR of AliFMDCorrector");
7e4038b5	79	TIter next(&o.fRingHistos);
	80	TObject* obj = 0;
	81	while ((obj = next())) fRingHistos.Add(obj);
	82	}
	83
	84	//____________________________________________________________________
72cc12cd	85	AliFMDCorrector::~AliFMDCorrector()
7e4038b5	86	{
7984e5f7	87	// Destructor
	88	//
	89	//
6ab100ec	90	DGUARD(fDebug, 3, "DTOR of AliFMDCorrector");
7e4038b5	91	fRingHistos.Delete();
	92	}
	93
	94	//____________________________________________________________________
72cc12cd	95	AliFMDCorrector&
72cc12cd	96	AliFMDCorrector::operator=(const AliFMDCorrector& o)
7e4038b5	97	{
7984e5f7	98	// Assignment operator
	99	//
	100	// Parameters:
	101	// o Object to assign from
6ab100ec	102	DGUARD(fDebug, 3, "Assignment of AliFMDCorrector");
d015ecfe	103	if (&o == this) return *this;
ea3e5d95	104	TNamed::operator=(o);
7e4038b5	105
ea3e5d95	106	fDebug = o.fDebug;
7e4038b5	107	fRingHistos.Delete();
7ec4d843	108	fUseSecondaryMap = o.fUseSecondaryMap;
	109	fUseVertexBias = o.fUseVertexBias;
	110	fUseAcceptance = o.fUseAcceptance;
81eda625	111	fUseMergingEfficiency = o.fUseMergingEfficiency;
7e4038b5	112	TIter next(&o.fRingHistos);
	113	TObject* obj = 0;
	114	while ((obj = next())) fRingHistos.Add(obj);
	115
	116	return *this;
	117	}
	118
7ec4d843	119	//____________________________________________________________________
7ec4d843	120	void
5934a3e3	121	AliFMDCorrector::SetupForData(const TAxis&)
7ec4d843	122	{
	123	//
	124	// Initialize this object
	125	//
	126	// Parameters:
	127	// etaAxis Eta axis to use
	128	//
6ab100ec	129	DGUARD(fDebug, 1, "Initialization of AliFMDCorrector");
7ec4d843	130	if (!fUseSecondaryMap)
	131	AliWarning("Secondary maps not used - BE CAREFUL");
	132	if (!fUseVertexBias)
	133	AliWarning("Vertex bias not used");
	134	if (!fUseAcceptance)
	135	AliWarning("Acceptance from dead-channels not used");
	136	}
	137
7e4038b5	138	//____________________________________________________________________
72cc12cd	139	AliFMDCorrector::RingHistos*
72cc12cd	140	AliFMDCorrector::GetRingHistos(UShort_t d, Char_t r) const
7e4038b5	141	{
7984e5f7	142	//
	143	// Get the ring histogram container
	144	// Parameters:
	145	// d Detector
	146	// r Ring
	147	//
	148	// Return:
	149	// Ring histogram container
	150	//
7e4038b5	151	Int_t idx = -1;
	152	switch (d) {
	153	case 1: idx = 0; break;
	154	case 2: idx = 1 + (r == 'I' \|\| r == 'i' ? 0 : 1); break;
	155	case 3: idx = 3 + (r == 'I' \|\| r == 'i' ? 0 : 1); break;
	156	}
	157	if (idx < 0 \|\| idx >= fRingHistos.GetEntries()) return 0;
	158
	159	return static_cast<RingHistos*>(fRingHistos.At(idx));
	160	}
	161
5ca83fee	162	//____________________________________________________________________
	163	void
	164	AliFMDCorrector::DivideMap(TH2* num, const TH2* denom,
	165	Bool_t alsoUnderOver) const
	166	{
	167	//
	168	// Implement TH1::Divide but
	169	// - Assume compatible histograms
	170	// - Unless third argument is true, do not divide over/under flow bins
	171	//
	172	if (!num \|\| !denom) return;
	173
	174	Int_t first = (alsoUnderOver ? 0 : 1);
	175	Int_t lastX = num->GetNbinsX() + (alsoUnderOver ? 1 : 0);
	176	Int_t lastY = num->GetNbinsY() + (alsoUnderOver ? 1 : 0);
	177
	178	for (Int_t ix = first; ix <= lastX; ix++) {
	179	for (Int_t iy = first; iy <= lastY; iy++) {
	180	Int_t bin = num->GetBin(ix,iy);
	181	Double_t c0 = num->GetBinContent(bin);
	182	Double_t c1 = denom->GetBinContent(bin);
	183	if (!c1) {
	184	num->SetBinContent(bin,0);
	185	num->SetBinError(bin, 0);
	186	continue;
	187	}
	188	Double_t w = c0 / c1;
	189	Double_t e0 = num->GetBinError(bin);
	190	Double_t e1 = denom->GetBinError(bin);
	191	Double_t c12 = c1*c1;
	192	Double_t e2 = (e0e0c1c1 + e1e1c0c0)/(c12*c12);
	193
	194	num->SetBinContent(bin, w);
	195	num->SetBinError(bin, TMath::Sqrt(e2));
	196	}
	197	}
	198	}
7e4038b5	199	//____________________________________________________________________
7e4038b5	200	Bool_t
72cc12cd	201	AliFMDCorrector::Correct(AliForwardUtil::Histos& hists,
7ec4d843	202	UShort_t vtxbin)
7e4038b5	203	{
7984e5f7	204	//
	205	// Do the calculations
	206	// Parameters:
	207	// hists Cache of histograms
	208	// vtxBin Vertex bin
	209	//
	210	// Return:
	211	// true on successs
	212	//
6ab100ec	213	DGUARD(fDebug, 1, "Correct histograms of AliFMDCorrector");
0bd4b00f	214	AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
7e4038b5	215
0bd4b00f	216	UShort_t uvb = vtxbin;
7e4038b5	217	for (UShort_t d=1; d<=3; d++) {
	218	UShort_t nr = (d == 1 ? 1 : 2);
	219	for (UShort_t q=0; q<nr; q++) {
72cc12cd	220	Char_t r = (q == 0 ? 'I' : 'O');
	221	TH2D* h = hists.Get(d,r);
	222	RingHistos* rh = GetRingHistos(d,r);
7ec4d843	223
	224	if (fUseSecondaryMap) {
	225	TH2D* bg = fcm.GetSecondaryMap()->GetCorrection(d, r, uvb);
	226	if (!bg) {
5bb5d1f6	227	AliWarning(Form("No secondary correction for FMDM%d%c "
5bb5d1f6	228	"in vertex bin %d", d, r, uvb));
7ec4d843	229	continue;
	230	}
	231	// Divide by primary/total ratio
5ca83fee	232	DivideMap(h, bg, false);
7e4038b5	233	}
7ec4d843	234	if (fUseVertexBias) {
	235	TH2D* ef = fcm.GetVertexBias()->GetCorrection(r, uvb);
	236	if (!ef) {
	237	AliWarning(Form("No event %s vertex bias correction in vertex bin %d",
	238	(r == 'I' \|\| r == 'i' ? "inner" : "outer"), uvb));
	239	continue;
	240	}
	241	// Divide by the event selection efficiency
5ca83fee	242	DivideMap(h, ef, false);
72cc12cd	243	}
7ec4d843	244	if (fUseAcceptance) {
	245	TH2D* ac = fcm.GetAcceptance()->GetCorrection(d, r, uvb);
	246	if (!ac) {
5bb5d1f6	247	AliWarning(Form("No acceptance correction for FMD%d%c in "
5bb5d1f6	248	"vertex bin %d", d, r, uvb));
7ec4d843	249	continue;
7ec4d843	250	}
5ca83fee	251	// Fill overflow bin with ones
	252	for (Int_t i = 1; i <= h->GetNbinsX(); i++)
	253	h->SetBinContent(i, h->GetNbinsY()+1, 1);
	254
	255	// Divide by the acceptance correction -
	256	DivideMap(h, ac, fcm.GetAcceptance()->HasOverflow());
7e4038b5	257	}
7e4038b5	258
81eda625	259	if (fUseMergingEfficiency) {
	260	if (!fcm.GetMergingEfficiency()) {
	261	AliWarning("No merging efficiencies");
	262	continue;
	263	}
	264	TH1D* sf = fcm.GetMergingEfficiency()->GetCorrection(d,r,uvb);
	265	if (!fcm.GetMergingEfficiency()->GetCorrection(d,r,uvb)) {
	266	AliWarning(Form("No merging efficiency for FMD%d%c at vertex bin %d",
	267	d, r, uvb));
	268	continue;
	269	}
7e4038b5	270
81eda625	271
	272	for (Int_t ieta = 1; ieta <= h->GetNbinsX(); ieta++) {
	273	Float_t c = sf->GetBinContent(ieta);
	274	Float_t ec = sf->GetBinError(ieta);
5ca83fee	275
81eda625	276	for (Int_t iphi = 1; iphi <= h->GetNbinsY(); iphi++) {
5ca83fee	277	if (c == 0) {
	278	h->SetBinContent(ieta,iphi,0);
	279	h->SetBinError(ieta,iphi,0);
	280	continue;
	281	}
	282
81eda625	283	Double_t m = h->GetBinContent(ieta, iphi) / c;
81eda625	284	Double_t em = h->GetBinError(ieta, iphi);
7e4038b5	285
81eda625	286	Double_t e = TMath::Sqrt(em * em + (m * ec) * (m * ec));
	287
	288	h->SetBinContent(ieta,iphi,m);
	289	h->SetBinError(ieta,iphi,e);
	290	}
7e4038b5	291	}
7e4038b5	292	}
12fffad7	293
7e4038b5	294	rh->fDensity->Add(h);
	295	}
	296	}
	297
	298	return kTRUE;
	299	}
	300
	301	//____________________________________________________________________
	302	void
5934a3e3	303	AliFMDCorrector::Terminate(const TList* dir, TList* output, Int_t nEvents)
7e4038b5	304	{
7984e5f7	305	//
	306	// Scale the histograms to the total number of events
	307	// Parameters:
	308	// dir Where the output is stored
	309	// nEvents Number of events
	310	//
6ab100ec	311	DGUARD(fDebug, 1, "Scale histograms of AliFMDCorrector");
7e4038b5	312	if (nEvents <= 0) return;
9d99b0dd	313	TList* d = static_cast<TList*>(dir->FindObject(GetName()));
9d99b0dd	314	if (!d) return;
7e4038b5	315
5934a3e3	316	TList* out = new TList;
	317	out->SetName(d->GetName());
	318	out->SetOwner();
	319
7e4038b5	320	TIter next(&fRingHistos);
7e4038b5	321	RingHistos* o = 0;
5bb5d1f6	322	THStack* sums = new THStack("sums", "Sums of ring results");
5bb5d1f6	323	while ((o = static_cast<RingHistos*>(next()))) {
5934a3e3	324	o->Terminate(d, nEvents);
	325	if (!o->fDensity) {
	326	Warning("Terminate", "No density from %s", o->GetName());
	327	continue;
	328	}
5bb5d1f6	329	TH1D* sum = o->fDensity->ProjectionX(o->GetName(), 1,
	330	o->fDensity->GetNbinsY(),"e");
	331	sum->Scale(1., "width");
	332	sum->SetTitle(o->GetName());
	333	sum->SetDirectory(0);
	334	sum->SetYTitle("#sum N_{ch,primary}");
	335	sums->Add(sum);
	336	}
5934a3e3	337	out->Add(sums);
5934a3e3	338	output->Add(out);
7e4038b5	339	}
7e4038b5	340	//____________________________________________________________________
7e4038b5	341	void
5934a3e3	342	AliFMDCorrector::CreateOutputObjects(TList* dir)
7e4038b5	343	{
fb3430ac	344	//
	345	// Output diagnostic histograms to directory
	346	//
	347	// Parameters:
	348	// dir List to write in
	349	//
6ab100ec	350	DGUARD(fDebug, 1, "Define output of AliFMDCorrector");
7e4038b5	351	TList* d = new TList;
	352	d->SetName(GetName());
	353	dir->Add(d);
241cca4d	354
	355	d->Add(AliForwardUtil::MakeParameter("secondary", fUseSecondaryMap));
	356	d->Add(AliForwardUtil::MakeParameter("vertexBias", fUseVertexBias));
	357	d->Add(AliForwardUtil::MakeParameter("acceptance", fUseAcceptance));
	358	d->Add(AliForwardUtil::MakeParameter("merging", fUseMergingEfficiency));
	359
	360
7e4038b5	361	TIter next(&fRingHistos);
	362	RingHistos* o = 0;
	363	while ((o = static_cast<RingHistos*>(next()))) {
5934a3e3	364	o->CreateOutputObjects(d);
7e4038b5	365	}
	366	}
	367
c8b1a7db	368	#define PF(N,V,...) \
	369	AliForwardUtil::PrintField(N,V, ## __VA_ARGS__)
	370	#define PFB(N,FLAG) \
	371	do { \
	372	AliForwardUtil::PrintName(N); \
	373	std::cout << std::boolalpha << (FLAG) << std::noboolalpha << std::endl; \
	374	} while(false)
	375	#define PFV(N,VALUE) \
	376	do { \
	377	AliForwardUtil::PrintName(N); \
	378	std::cout << (VALUE) << std::endl; } while(false)
	379
0bd4b00f	380	//____________________________________________________________________
0bd4b00f	381	void
72cc12cd	382	AliFMDCorrector::Print(Option_t* /* option */) const
0bd4b00f	383	{
7984e5f7	384	//
	385	// Print information
	386	// Parameters:
	387	// option Not used
	388	//
c8b1a7db	389	AliForwardUtil::PrintTask(*this);
	390	gROOT->IncreaseDirLevel();
	391	PFB("Use secondary maps", fUseSecondaryMap );
	392	PFB("Use vertex bias", fUseVertexBias );
	393	PFB("Use acceptance", fUseAcceptance );
	394	PFB("Use merging efficiency", fUseMergingEfficiency);
	395	gROOT->DecreaseDirLevel();
0bd4b00f	396	}
0bd4b00f	397
7e4038b5	398	//====================================================================
72cc12cd	399	AliFMDCorrector::RingHistos::RingHistos()
9d99b0dd	400	: AliForwardUtil::RingHistos(),
7e4038b5	401	fDensity(0)
7984e5f7	402	{
	403	// Constructor
	404	//
	405	//
	406	}
7e4038b5	407
7e4038b5	408	//____________________________________________________________________
72cc12cd	409	AliFMDCorrector::RingHistos::RingHistos(UShort_t d, Char_t r)
9d99b0dd	410	: AliForwardUtil::RingHistos(d,r),
7e4038b5	411	fDensity(0)
7e4038b5	412	{
7984e5f7	413	//
	414	// Constructor
	415	// Parameters:
	416	// d detector
	417	// r ring
	418	//
5bb5d1f6	419	fDensity = new TH2D("primaryDensity",
5bb5d1f6	420	"#sum N_{ch,primary}/(#Delta#eta#Delta#phi)",
7e4038b5	421	200, -4, 6, (r == 'I' \|\| r == 'i' ? 20 : 40),
	422	0, 2*TMath::Pi());
	423	fDensity->SetDirectory(0);
5bb5d1f6	424	fDensity->SetMarkerColor(Color());
5bb5d1f6	425	fDensity->Sumw2();
7e4038b5	426	fDensity->SetXTitle("#eta");
	427	fDensity->SetYTitle("#phi [radians]");
	428	fDensity->SetZTitle("Primary N_{ch} density");
	429	}
	430	//____________________________________________________________________
72cc12cd	431	AliFMDCorrector::RingHistos::RingHistos(const RingHistos& o)
9d99b0dd	432	: AliForwardUtil::RingHistos(o),
7e4038b5	433	fDensity(o.fDensity)
7984e5f7	434	{
	435	//
	436	// Copy constructor
	437	// Parameters:
	438	// o Object to copy from
	439	//
	440	}
7e4038b5	441
7e4038b5	442	//____________________________________________________________________
72cc12cd	443	AliFMDCorrector::RingHistos&
72cc12cd	444	AliFMDCorrector::RingHistos::operator=(const RingHistos& o)
7e4038b5	445	{
7984e5f7	446	//
	447	// Assignment operator
	448	// Parameters:
	449	// o Object to assign from
	450	//
	451	// Return:
	452	// Reference to this
	453	//
d015ecfe	454	if (&o == this) return *this;
9d99b0dd	455	AliForwardUtil::RingHistos::operator=(o);
7e4038b5	456
	457	if (fDensity) delete fDensity;
	458
	459	fDensity = static_cast<TH2D*>(o.fDensity->Clone());
	460
	461	return *this;
	462	}
	463	//____________________________________________________________________
72cc12cd	464	AliFMDCorrector::RingHistos::~RingHistos()
7e4038b5	465	{
7984e5f7	466	//
	467	// Destructor
	468	//
b7ab8a2c	469	// if (fDensity) delete fDensity;
7e4038b5	470	}
	471
	472	//____________________________________________________________________
	473	void
5934a3e3	474	AliFMDCorrector::RingHistos::CreateOutputObjects(TList* dir)
7e4038b5	475	{
7984e5f7	476	//
	477	// Make output
	478	// Parameters:
	479	// dir Where to put it
	480	//
9d99b0dd	481	TList* d = DefineOutputList(dir);
7e4038b5	482	d->Add(fDensity);
9d99b0dd	483	}
	484
	485	//____________________________________________________________________
	486	void
5934a3e3	487	AliFMDCorrector::RingHistos::Terminate(TList* dir, Int_t nEvents)
9d99b0dd	488	{
7984e5f7	489	//
	490	// Scale the histograms to the total number of events
	491	// Parameters:
	492	// dir where the output is stored
	493	// nEvents Number of events
	494	//
9d99b0dd	495	TList* l = GetOutputList(dir);
	496	if (!l) return;
	497
5934a3e3	498	TH2D* density = static_cast<TH2D*>(GetOutputHist(l,"primaryDensity"));
5bb5d1f6	499	if (density) density->Scale(1./nEvents);
5934a3e3	500	fDensity = density;
7e4038b5	501	}
	502
	503	//____________________________________________________________________
	504	//
	505	// EOF
	506	//
	507
	508
	509