[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 "AliFMDCorrDoubleHit.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
}

//____________________________________________________________________
AliFMDCorrector::AliFMDCorrector(const char* title)
  : TNamed("fmdCorrector", title), 
    fRingHistos(), 
    fUseSecondaryMap(true),
    fUseVertexBias(false),
    fUseAcceptance(true),
    fUseMergingEfficiency(false),
    fDebug(0)
{
  // Constructor 
  // 
  // Parameters: 
  //   title   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 
  TIter    next(&o.fRingHistos);
  TObject* obj = 0;
  while ((obj = next())) fRingHistos.Add(obj);
}

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

//____________________________________________________________________
AliFMDCorrector&
AliFMDCorrector::operator=(const AliFMDCorrector& o)
{
  // Assignment operator 
  // 
  // Parameters:
  //   o   Object to assign from 
  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::Init(const TAxis&)
{
  //
  // Initialize this object
  //
  // Parameters:
  //    etaAxis Eta axis to use
  //
  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));
}
    
//____________________________________________________________________
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 
  //
  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
        h->Divide(bg);
      }
      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
        h->Divide(ef);
      }
      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;
        }
        // Divide by the acceptance correction
        h->Divide(ac);
      }

      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);
	  
	  if (c == 0) continue;
	  
	  for (Int_t iphi = 1; iphi <= h->GetNbinsY(); iphi++) { 
	    Double_t m  = h->GetBinContent(ieta, iphi) / c;
	    Double_t em = h->GetBinError(ieta, iphi);
	  
	    Double_t e  = TMath::Sqrt(em * em + (m * ec) * (m * ec));
	    
	    h->SetBinContent(ieta,iphi,m);
	    h->SetBinError(ieta,iphi,e);
	  }
	}
      }
      //HHD
      /*
      TH2D*  bg = fcm.GetSecondaryMap()->GetCorrection(d, r, uvb);
      TH2D   hRing("hring","hring",bg->GetNbinsX(),
		   bg->GetXaxis()->GetXmin(),
		   bg->GetXaxis()->GetXmax(),
		   bg->GetNbinsY(),
		   bg->GetYaxis()->GetXmin(),
		   bg->GetYaxis()->GetXmax());
      
      Int_t edgebin[4] = {0,0,0,0};
      for(Int_t ii = 1; ii <=bg->GetNbinsX(); ii++) {
	for(Int_t jj = 1; jj <=bg->GetNbinsY(); jj++) {
	  Float_t bgcor = bg->GetBinContent(ii,jj);
	  if(bgcor<0.1) continue;
	  if(edgebin[0] == 0) edgebin[0] = ii;
	  if(edgebin[0] == ii) continue;
	  if(edgebin[0] > 0 && edgebin[1] == 0) edgebin[1] = ii;
	  if(edgebin[0]>0 && edgebin[1]>0) break; 
	}
      }
      for(Int_t ii = bg->GetNbinsX(); ii >= 1;  ii--) {
	for(Int_t jj = 1; jj <=bg->GetNbinsY(); jj++) {
	  Float_t bgcor = bg->GetBinContent(ii,jj);
	  if(bgcor<0.1) continue;
	  if(edgebin[2] == 0) edgebin[2] = ii;
	  if(edgebin[2] == ii) continue;
	  if(edgebin[2] > 0 && edgebin[3] == 0) edgebin[3] = ii;
	  if(edgebin[2]>0 && edgebin[3]>0) break; 
	}
      }
      for(Int_t ii = 1; ii <=bg->GetNbinsX(); ii++) {
	for(Int_t jj = 1; jj <=bg->GetNbinsY(); jj++) {
	  Float_t data = h->GetBinContent(ii,jj);
	  if(data <0.000001) continue;
	  if(edgebin[0] == ii || edgebin[1] == ii || edgebin[2] == ii || edgebin[3] == ii) continue;
	  hRing.SetBinContent(ii,jj,data);
	  hRing.SetBinError(ii,jj,h->GetBinError(ii,jj));
	}
      }
      
      //std::cout<<edgebin[0]<<"  "<<edgebin[1]<<"  "<<edgebin[2]<<"   "<<edgebin[3]<<std::endl;
      */
      
      rh->fDensity->Add(h);
    }
  }
  
  return kTRUE;
}

//____________________________________________________________________
void
AliFMDCorrector::ScaleHistograms(const 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 
  //
  if (nEvents <= 0) return;
  TList* d = static_cast<TList*>(dir->FindObject(GetName()));
  if (!d) return;

  TIter    next(&fRingHistos);
  RingHistos* o = 0;
  THStack* sums = new THStack("sums", "Sums of ring results");
  while ((o = static_cast<RingHistos*>(next()))) {
    o->ScaleHistograms(d, nEvents);
    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);
  }
  d->Add(sums);

}
//____________________________________________________________________
void
AliFMDCorrector::DefineOutput(TList* dir)
{
  // 
  // Output diagnostic histograms to directory 
  // 
  // Parameters:
  //    dir List to write in
  //  
  TList* d = new TList;
  d->SetName(GetName());
  dir->Add(d);
  TIter    next(&fRingHistos);
  RingHistos* o = 0;
  while ((o = static_cast<RingHistos*>(next()))) {
    o->Output(d);
  }
}

//____________________________________________________________________
void
AliFMDCorrector::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 << ClassName() << ": " << GetName() <<  "\n"
            << std::boolalpha
            << ind << " Use secondary maps:     " << fUseSecondaryMap << "\n"
            << ind << " Use vertex bias:        " << fUseVertexBias << "\n"
            << ind << " Use acceptance:         " << fUseAcceptance << "\n"
            << ind << " Use merging efficiency: " << fUseMergingEfficiency
            << std::noboolalpha << std::endl;
}

//====================================================================
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::Output(TList* dir)
{
  // 
  // Make output 
  // Parameters:
  //    dir Where to put it 
  //
  TList* d = DefineOutputList(dir);
  d->Add(fDensity);
}

//____________________________________________________________________
void
AliFMDCorrector::RingHistos::ScaleHistograms(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; 

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

//____________________________________________________________________
//
// 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"
fb3430ac	11	// #include "AliFMDCorrDoubleHit.h"
	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	{
	37	// Constructor
	38	}
7e4038b5	39
7e4038b5	40	//____________________________________________________________________
72cc12cd	41	AliFMDCorrector::AliFMDCorrector(const char* title)
cc83fca2	42	: TNamed("fmdCorrector", title),
7e4038b5	43	fRingHistos(),
7ec4d843	44	fUseSecondaryMap(true),
9d05ffeb	45	fUseVertexBias(false),
7ec4d843	46	fUseAcceptance(true),
9d05ffeb	47	fUseMergingEfficiency(false),
ea3e5d95	48	fDebug(0)
7e4038b5	49	{
7984e5f7	50	// Constructor
	51	//
	52	// Parameters:
	53	// title Title
7e4038b5	54	fRingHistos.SetName(GetName());
	55	fRingHistos.Add(new RingHistos(1, 'I'));
	56	fRingHistos.Add(new RingHistos(2, 'I'));
	57	fRingHistos.Add(new RingHistos(2, 'O'));
	58	fRingHistos.Add(new RingHistos(3, 'I'));
	59	fRingHistos.Add(new RingHistos(3, 'O'));
	60	}
	61
	62	//____________________________________________________________________
72cc12cd	63	AliFMDCorrector::AliFMDCorrector(const AliFMDCorrector& o)
7e4038b5	64	: TNamed(o),
7e4038b5	65	fRingHistos(),
7ec4d843	66	fUseSecondaryMap(o.fUseSecondaryMap),
	67	fUseVertexBias(o.fUseVertexBias),
	68	fUseAcceptance(o.fUseAcceptance),
81eda625	69	fUseMergingEfficiency(o.fUseMergingEfficiency),
ea3e5d95	70	fDebug(o.fDebug)
7e4038b5	71	{
7984e5f7	72	// Copy constructor
	73	//
	74	// Parameters:
	75	// o Object to copy from
7e4038b5	76	TIter next(&o.fRingHistos);
	77	TObject* obj = 0;
	78	while ((obj = next())) fRingHistos.Add(obj);
	79	}
	80
	81	//____________________________________________________________________
72cc12cd	82	AliFMDCorrector::~AliFMDCorrector()
7e4038b5	83	{
7984e5f7	84	// Destructor
	85	//
	86	//
7e4038b5	87	fRingHistos.Delete();
	88	}
	89
	90	//____________________________________________________________________
72cc12cd	91	AliFMDCorrector&
72cc12cd	92	AliFMDCorrector::operator=(const AliFMDCorrector& o)
7e4038b5	93	{
7984e5f7	94	// Assignment operator
	95	//
	96	// Parameters:
	97	// o Object to assign from
d015ecfe	98	if (&o == this) return *this;
ea3e5d95	99	TNamed::operator=(o);
7e4038b5	100
ea3e5d95	101	fDebug = o.fDebug;
7e4038b5	102	fRingHistos.Delete();
7ec4d843	103	fUseSecondaryMap = o.fUseSecondaryMap;
	104	fUseVertexBias = o.fUseVertexBias;
	105	fUseAcceptance = o.fUseAcceptance;
81eda625	106	fUseMergingEfficiency = o.fUseMergingEfficiency;
7e4038b5	107	TIter next(&o.fRingHistos);
	108	TObject* obj = 0;
	109	while ((obj = next())) fRingHistos.Add(obj);
	110
	111	return *this;
	112	}
	113
7ec4d843	114	//____________________________________________________________________
	115	void
	116	AliFMDCorrector::Init(const TAxis&)
	117	{
	118	//
	119	// Initialize this object
	120	//
	121	// Parameters:
	122	// etaAxis Eta axis to use
	123	//
	124	if (!fUseSecondaryMap)
	125	AliWarning("Secondary maps not used - BE CAREFUL");
	126	if (!fUseVertexBias)
	127	AliWarning("Vertex bias not used");
	128	if (!fUseAcceptance)
	129	AliWarning("Acceptance from dead-channels not used");
	130	}
	131
7e4038b5	132	//____________________________________________________________________
72cc12cd	133	AliFMDCorrector::RingHistos*
72cc12cd	134	AliFMDCorrector::GetRingHistos(UShort_t d, Char_t r) const
7e4038b5	135	{
7984e5f7	136	//
	137	// Get the ring histogram container
	138	// Parameters:
	139	// d Detector
	140	// r Ring
	141	//
	142	// Return:
	143	// Ring histogram container
	144	//
7e4038b5	145	Int_t idx = -1;
	146	switch (d) {
	147	case 1: idx = 0; break;
	148	case 2: idx = 1 + (r == 'I' \|\| r == 'i' ? 0 : 1); break;
	149	case 3: idx = 3 + (r == 'I' \|\| r == 'i' ? 0 : 1); break;
	150	}
	151	if (idx < 0 \|\| idx >= fRingHistos.GetEntries()) return 0;
	152
	153	return static_cast<RingHistos*>(fRingHistos.At(idx));
	154	}
	155
	156	//____________________________________________________________________
	157	Bool_t
72cc12cd	158	AliFMDCorrector::Correct(AliForwardUtil::Histos& hists,
7ec4d843	159	UShort_t vtxbin)
7e4038b5	160	{
7984e5f7	161	//
	162	// Do the calculations
	163	// Parameters:
	164	// hists Cache of histograms
	165	// vtxBin Vertex bin
	166	//
	167	// Return:
	168	// true on successs
	169	//
0bd4b00f	170	AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
7e4038b5	171
0bd4b00f	172	UShort_t uvb = vtxbin;
7e4038b5	173	for (UShort_t d=1; d<=3; d++) {
	174	UShort_t nr = (d == 1 ? 1 : 2);
	175	for (UShort_t q=0; q<nr; q++) {
72cc12cd	176	Char_t r = (q == 0 ? 'I' : 'O');
	177	TH2D* h = hists.Get(d,r);
	178	RingHistos* rh = GetRingHistos(d,r);
7ec4d843	179
	180	if (fUseSecondaryMap) {
	181	TH2D* bg = fcm.GetSecondaryMap()->GetCorrection(d, r, uvb);
	182	if (!bg) {
5bb5d1f6	183	AliWarning(Form("No secondary correction for FMDM%d%c "
5bb5d1f6	184	"in vertex bin %d", d, r, uvb));
7ec4d843	185	continue;
	186	}
	187	// Divide by primary/total ratio
	188	h->Divide(bg);
7e4038b5	189	}
7ec4d843	190	if (fUseVertexBias) {
	191	TH2D* ef = fcm.GetVertexBias()->GetCorrection(r, uvb);
	192	if (!ef) {
	193	AliWarning(Form("No event %s vertex bias correction in vertex bin %d",
	194	(r == 'I' \|\| r == 'i' ? "inner" : "outer"), uvb));
	195	continue;
	196	}
	197	// Divide by the event selection efficiency
	198	h->Divide(ef);
72cc12cd	199	}
7ec4d843	200	if (fUseAcceptance) {
	201	TH2D* ac = fcm.GetAcceptance()->GetCorrection(d, r, uvb);
	202	if (!ac) {
5bb5d1f6	203	AliWarning(Form("No acceptance correction for FMD%d%c in "
5bb5d1f6	204	"vertex bin %d", d, r, uvb));
7ec4d843	205	continue;
	206	}
	207	// Divide by the acceptance correction
	208	h->Divide(ac);
7e4038b5	209	}
7e4038b5	210
81eda625	211	if (fUseMergingEfficiency) {
	212	if (!fcm.GetMergingEfficiency()) {
	213	AliWarning("No merging efficiencies");
	214	continue;
	215	}
	216	TH1D* sf = fcm.GetMergingEfficiency()->GetCorrection(d,r,uvb);
	217	if (!fcm.GetMergingEfficiency()->GetCorrection(d,r,uvb)) {
	218	AliWarning(Form("No merging efficiency for FMD%d%c at vertex bin %d",
	219	d, r, uvb));
	220	continue;
	221	}
7e4038b5	222
81eda625	223
	224	for (Int_t ieta = 1; ieta <= h->GetNbinsX(); ieta++) {
	225	Float_t c = sf->GetBinContent(ieta);
	226	Float_t ec = sf->GetBinError(ieta);
	227
	228	if (c == 0) continue;
7e4038b5	229
81eda625	230	for (Int_t iphi = 1; iphi <= h->GetNbinsY(); iphi++) {
	231	Double_t m = h->GetBinContent(ieta, iphi) / c;
	232	Double_t em = h->GetBinError(ieta, iphi);
7e4038b5	233
81eda625	234	Double_t e = TMath::Sqrt(em * em + (m * ec) * (m * ec));
	235
	236	h->SetBinContent(ieta,iphi,m);
	237	h->SetBinError(ieta,iphi,e);
	238	}
7e4038b5	239	}
7e4038b5	240	}
12fffad7	241	//HHD
	242	/*
	243	TH2D* bg = fcm.GetSecondaryMap()->GetCorrection(d, r, uvb);
	244	TH2D hRing("hring","hring",bg->GetNbinsX(),
	245	bg->GetXaxis()->GetXmin(),
	246	bg->GetXaxis()->GetXmax(),
	247	bg->GetNbinsY(),
	248	bg->GetYaxis()->GetXmin(),
	249	bg->GetYaxis()->GetXmax());
	250
	251	Int_t edgebin[4] = {0,0,0,0};
	252	for(Int_t ii = 1; ii <=bg->GetNbinsX(); ii++) {
	253	for(Int_t jj = 1; jj <=bg->GetNbinsY(); jj++) {
	254	Float_t bgcor = bg->GetBinContent(ii,jj);
	255	if(bgcor<0.1) continue;
	256	if(edgebin[0] == 0) edgebin[0] = ii;
	257	if(edgebin[0] == ii) continue;
	258	if(edgebin[0] > 0 && edgebin[1] == 0) edgebin[1] = ii;
	259	if(edgebin[0]>0 && edgebin[1]>0) break;
	260	}
	261	}
	262	for(Int_t ii = bg->GetNbinsX(); ii >= 1; ii--) {
	263	for(Int_t jj = 1; jj <=bg->GetNbinsY(); jj++) {
	264	Float_t bgcor = bg->GetBinContent(ii,jj);
	265	if(bgcor<0.1) continue;
	266	if(edgebin[2] == 0) edgebin[2] = ii;
	267	if(edgebin[2] == ii) continue;
	268	if(edgebin[2] > 0 && edgebin[3] == 0) edgebin[3] = ii;
	269	if(edgebin[2]>0 && edgebin[3]>0) break;
	270	}
	271	}
	272	for(Int_t ii = 1; ii <=bg->GetNbinsX(); ii++) {
	273	for(Int_t jj = 1; jj <=bg->GetNbinsY(); jj++) {
	274	Float_t data = h->GetBinContent(ii,jj);
	275	if(data <0.000001) continue;
	276	if(edgebin[0] == ii \|\| edgebin[1] == ii \|\| edgebin[2] == ii \|\| edgebin[3] == ii) continue;
	277	hRing.SetBinContent(ii,jj,data);
	278	hRing.SetBinError(ii,jj,h->GetBinError(ii,jj));
	279	}
	280	}
	281
	282	//std::cout<<edgebin[0]<<" "<<edgebin[1]<<" "<<edgebin[2]<<" "<<edgebin[3]<<std::endl;
	283	*/
	284
7e4038b5	285	rh->fDensity->Add(h);
	286	}
	287	}
	288
	289	return kTRUE;
	290	}
	291
	292	//____________________________________________________________________
	293	void
fb3430ac	294	AliFMDCorrector::ScaleHistograms(const TList* dir, Int_t nEvents)
7e4038b5	295	{
7984e5f7	296	//
	297	// Scale the histograms to the total number of events
	298	// Parameters:
	299	// dir Where the output is stored
	300	// nEvents Number of events
	301	//
7e4038b5	302	if (nEvents <= 0) return;
9d99b0dd	303	TList* d = static_cast<TList*>(dir->FindObject(GetName()));
9d99b0dd	304	if (!d) return;
7e4038b5	305
	306	TIter next(&fRingHistos);
	307	RingHistos* o = 0;
5bb5d1f6	308	THStack* sums = new THStack("sums", "Sums of ring results");
5bb5d1f6	309	while ((o = static_cast<RingHistos*>(next()))) {
9d99b0dd	310	o->ScaleHistograms(d, nEvents);
5bb5d1f6	311	TH1D* sum = o->fDensity->ProjectionX(o->GetName(), 1,
	312	o->fDensity->GetNbinsY(),"e");
	313	sum->Scale(1., "width");
	314	sum->SetTitle(o->GetName());
	315	sum->SetDirectory(0);
	316	sum->SetYTitle("#sum N_{ch,primary}");
	317	sums->Add(sum);
	318	}
	319	d->Add(sums);
	320
7e4038b5	321	}
7e4038b5	322	//____________________________________________________________________
7e4038b5	323	void
72cc12cd	324	AliFMDCorrector::DefineOutput(TList* dir)
7e4038b5	325	{
fb3430ac	326	//
	327	// Output diagnostic histograms to directory
	328	//
	329	// Parameters:
	330	// dir List to write in
	331	//
7e4038b5	332	TList* d = new TList;
	333	d->SetName(GetName());
	334	dir->Add(d);
	335	TIter next(&fRingHistos);
	336	RingHistos* o = 0;
	337	while ((o = static_cast<RingHistos*>(next()))) {
	338	o->Output(d);
	339	}
	340	}
	341
0bd4b00f	342	//____________________________________________________________________
0bd4b00f	343	void
72cc12cd	344	AliFMDCorrector::Print(Option_t* /* option */) const
0bd4b00f	345	{
7984e5f7	346	//
	347	// Print information
	348	// Parameters:
	349	// option Not used
	350	//
0bd4b00f	351	char ind[gROOT->GetDirLevel()+1];
	352	for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
	353	ind[gROOT->GetDirLevel()] = '\0';
1f480471	354	std::cout << ind << ClassName() << ": " << GetName() << "\n"
7ec4d843	355	<< std::boolalpha
	356	<< ind << " Use secondary maps: " << fUseSecondaryMap << "\n"
	357	<< ind << " Use vertex bias: " << fUseVertexBias << "\n"
	358	<< ind << " Use acceptance: " << fUseAcceptance << "\n"
	359	<< ind << " Use merging efficiency: " << fUseMergingEfficiency
8e2bb72a	360	<< std::noboolalpha << std::endl;
0bd4b00f	361	}
0bd4b00f	362
7e4038b5	363	//====================================================================
72cc12cd	364	AliFMDCorrector::RingHistos::RingHistos()
9d99b0dd	365	: AliForwardUtil::RingHistos(),
7e4038b5	366	fDensity(0)
7984e5f7	367	{
	368	// Constructor
	369	//
	370	//
	371	}
7e4038b5	372
7e4038b5	373	//____________________________________________________________________
72cc12cd	374	AliFMDCorrector::RingHistos::RingHistos(UShort_t d, Char_t r)
9d99b0dd	375	: AliForwardUtil::RingHistos(d,r),
7e4038b5	376	fDensity(0)
7e4038b5	377	{
7984e5f7	378	//
	379	// Constructor
	380	// Parameters:
	381	// d detector
	382	// r ring
	383	//
5bb5d1f6	384	fDensity = new TH2D("primaryDensity",
5bb5d1f6	385	"#sum N_{ch,primary}/(#Delta#eta#Delta#phi)",
7e4038b5	386	200, -4, 6, (r == 'I' \|\| r == 'i' ? 20 : 40),
	387	0, 2*TMath::Pi());
	388	fDensity->SetDirectory(0);
5bb5d1f6	389	fDensity->SetMarkerColor(Color());
5bb5d1f6	390	fDensity->Sumw2();
7e4038b5	391	fDensity->SetXTitle("#eta");
	392	fDensity->SetYTitle("#phi [radians]");
	393	fDensity->SetZTitle("Primary N_{ch} density");
	394	}
	395	//____________________________________________________________________
72cc12cd	396	AliFMDCorrector::RingHistos::RingHistos(const RingHistos& o)
9d99b0dd	397	: AliForwardUtil::RingHistos(o),
7e4038b5	398	fDensity(o.fDensity)
7984e5f7	399	{
	400	//
	401	// Copy constructor
	402	// Parameters:
	403	// o Object to copy from
	404	//
	405	}
7e4038b5	406
7e4038b5	407	//____________________________________________________________________
72cc12cd	408	AliFMDCorrector::RingHistos&
72cc12cd	409	AliFMDCorrector::RingHistos::operator=(const RingHistos& o)
7e4038b5	410	{
7984e5f7	411	//
	412	// Assignment operator
	413	// Parameters:
	414	// o Object to assign from
	415	//
	416	// Return:
	417	// Reference to this
	418	//
d015ecfe	419	if (&o == this) return *this;
9d99b0dd	420	AliForwardUtil::RingHistos::operator=(o);
7e4038b5	421
	422	if (fDensity) delete fDensity;
	423
	424	fDensity = static_cast<TH2D*>(o.fDensity->Clone());
	425
	426	return *this;
	427	}
	428	//____________________________________________________________________
72cc12cd	429	AliFMDCorrector::RingHistos::~RingHistos()
7e4038b5	430	{
7984e5f7	431	//
	432	// Destructor
	433	//
7e4038b5	434	if (fDensity) delete fDensity;
	435	}
	436
	437	//____________________________________________________________________
	438	void
72cc12cd	439	AliFMDCorrector::RingHistos::Output(TList* dir)
7e4038b5	440	{
7984e5f7	441	//
	442	// Make output
	443	// Parameters:
	444	// dir Where to put it
	445	//
9d99b0dd	446	TList* d = DefineOutputList(dir);
7e4038b5	447	d->Add(fDensity);
9d99b0dd	448	}
	449
	450	//____________________________________________________________________
	451	void
5bb5d1f6	452	AliFMDCorrector::RingHistos::ScaleHistograms(TList* dir, Int_t nEvents)
9d99b0dd	453	{
7984e5f7	454	//
	455	// Scale the histograms to the total number of events
	456	// Parameters:
	457	// dir where the output is stored
	458	// nEvents Number of events
	459	//
9d99b0dd	460	TList* l = GetOutputList(dir);
	461	if (!l) return;
	462
5bb5d1f6	463	TH1* density = GetOutputHist(l,"primaryDensity");
5bb5d1f6	464	if (density) density->Scale(1./nEvents);
7e4038b5	465	}
	466
	467	//____________________________________________________________________
	468	//
	469	// EOF
	470	//
	471
	472
	473