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

// 
// This class calculates the exclusive charged particle density
// in each for the 5 FMD rings. 
//
// Input:
//   - 5 RingHistos objects - each with a number of vertex dependent 
//     2D histograms of the inclusive charge particle density 
//
// Output:
//   - 5 RingHistos objects - each with a number of vertex dependent 
//     2D histograms of the exclusive charge particle density 
// 
// Corrections used: 
//   - AliFMDCorrSecondaryMap;
//   - AliFMDCorrVertexBias
//   - AliFMDCorrMergingEfficiency
//
#include "AliFMDMCCorrector.h"
#include <AliESDFMD.h>
#include <TAxis.h>
#include <TList.h>
#include <TMath.h>
#include "AliForwardCorrectionManager.h"
#include "AliFMDCorrVertexBias.h"
#include "AliLog.h"
#include <TH2D.h>
#include <TROOT.h>
#include <TProfile2D.h>
#include <iostream>
ClassImp(AliFMDMCCorrector)
#if 0
; // For Emacs
#endif 


//____________________________________________________________________
AliFMDMCCorrector::~AliFMDMCCorrector()
{
  // 
  // Destructor 
  //
  if (fComps) fComps->Clear();
  if (fFMD1i) delete fFMD1i;
  if (fFMD2i) delete fFMD2i;
  if (fFMD2o) delete fFMD2o;
  if (fFMD3i) delete fFMD3i;
  if (fFMD3o) delete fFMD3o;
}

//____________________________________________________________________
AliFMDMCCorrector&
AliFMDMCCorrector::operator=(const AliFMDMCCorrector& o)
{
  // 
  // Assignement operator
  // 
  // Parameters:
  //    o Object to assign from 
  // 
  // Return:
  //    Reference to this object
  //
  if (&o == this) return *this; 
  AliFMDCorrector::operator=(o);
  fSecondaryForMC = o.fSecondaryForMC;
  return *this;
}

//____________________________________________________________________
Bool_t
AliFMDMCCorrector::CorrectMC(AliForwardUtil::Histos& hists,
			       UShort_t                vtxbin)
{
  // 
  // Do the calculations 
  // 
  // Parameters:
  //    hists    Cache of histograms 
  //    vtxBin   Vertex bin 
  // 
  // Return:
  //    true on successs 
  //
  if ((!fUseSecondaryMap || !fSecondaryForMC) && fUseVertexBias) 
    return kTRUE;

  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);


      if (fUseSecondaryMap && fSecondaryForMC) {
        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 vertex bias correction in vertex bin %d",
			uvb));
          continue;
        }
        // Divide by the event selection efficiency
        h->Divide(ef);
      }
    }
  }
  
  return kTRUE;
}

//____________________________________________________________________
void
AliFMDMCCorrector::Init(const TAxis& eAxis)
{
  // 
  // Initialize this object 
  // 
  // Parameters:
  //    etaAxis Eta axis to use 
  //
  AliFMDCorrector::Init(eAxis);

  fFMD1i = Make(1,'I',eAxis);
  fFMD2i = Make(2,'I',eAxis);
  fFMD2o = Make(2,'O',eAxis);
  fFMD3i = Make(3,'I',eAxis);
  fFMD3o = Make(3,'O',eAxis);

  fComps->Add(fFMD1i);
  fComps->Add(fFMD2i);
  fComps->Add(fFMD2o);
  fComps->Add(fFMD3i);
  fComps->Add(fFMD3o);
}

//____________________________________________________________________
TProfile2D*
AliFMDMCCorrector::Make(UShort_t d, Char_t r, 
				const TAxis& axis) const
{
  // 
  // MAke comparison profiles
  // 
  // Parameters:
  //    d     Detector 
  //    r     Ring 
  //    axis  Eta axis 
  // 
  // Return:
  //    Newly allocated profile object
  //
  TProfile2D* ret = new TProfile2D(Form("FMD%d%c_esd_vs_mc", d, r),
				   Form("ESD/MC signal for FMD%d%c", d, r),
				   axis.GetNbins(), 
				   axis.GetXmin(),
				   axis.GetXmax(), 
				   (r == 'I' || r == 'i') ? 20 : 40,
				   0, 2*TMath::Pi());
  ret->GetXaxis()->SetTitle("#eta");
  ret->GetYaxis()->SetTitle("#varphi [degrees]");
  ret->GetZaxis()->SetTitle("#LT primary density ESD/MC#GT");
  ret->SetDirectory(0);
  return ret;
}
//____________________________________________________________________
void
AliFMDMCCorrector::Fill(UShort_t d, Char_t r, TH2* esd, TH2* mc)
{
  // 
  // Fill comparison profiles
  // 
  // Parameters:
  //    d    Detector 
  //    r    Ring 
  //    esd  ESD histogram
  //    mc   MC histogram
  //
  if (!esd || !mc) return;
  TProfile2D* p = 0;
  switch (d) { 
  case 1:  p = fFMD1i;                                   break;
  case 2:  p = (r == 'I' || r == 'i' ? fFMD2i : fFMD2o); break;
  case 3:  p = (r == 'I' || r == 'i' ? fFMD3i : fFMD3o); break;
  }
  if (!p) return;

  for (Int_t iEta = 1; iEta <= esd->GetNbinsX(); iEta++) { 
    Double_t eta = esd->GetXaxis()->GetBinCenter(iEta);
    for (Int_t iPhi = 1; iPhi <= esd->GetNbinsY(); iPhi++) { 
      Double_t phi  = esd->GetYaxis()->GetBinCenter(iPhi);
      Double_t mEsd = esd->GetBinContent(iEta,iPhi);
      Double_t mMc  = mc->GetBinContent(iEta,iPhi);
      
      p->Fill(eta, phi, (mMc > 0 ? mEsd / mMc : 0));
    }
  }
}

//____________________________________________________________________
Bool_t
AliFMDMCCorrector::CompareResults(AliForwardUtil::Histos& esd,
					  AliForwardUtil::Histos& mc)
{
  // 
  // Compare the result of analysing the ESD for 
  // the inclusive charged particle density to analysing 
  // MC truth 
  // 
  // Parameters:
  //    esd 
  //    mc 
  // 
  // Return:
  //   true 
  //
  Fill(1, 'I', esd.Get(1,'I'), mc.Get(1,'I'));
  Fill(2, 'I', esd.Get(2,'I'), mc.Get(2,'I'));
  Fill(2, 'O', esd.Get(2,'O'), mc.Get(2,'O'));
  Fill(3, 'I', esd.Get(3,'I'), mc.Get(3,'I'));
  Fill(3, 'O', esd.Get(3,'O'), mc.Get(3,'O'));

  return kTRUE;
}

//____________________________________________________________________
void
AliFMDMCCorrector::DefineOutput(TList* dir)
{
  // 
  // Output diagnostic histograms to directory 
  // 
  // Parameters:
  //    dir List to write in
  //  
  AliFMDCorrector::DefineOutput(dir);
  TList* d = static_cast<TList*>(dir->FindObject(GetName()));

  fComps = new TList;
  fComps->SetName("esd_mc_comparison");
  d->Add(fComps);
}
//____________________________________________________________________
void
AliFMDMCCorrector::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';
  AliFMDCorrector::Print(option);
  std::cout << std::boolalpha
            << ind << " Use sec. map on MC:     " << fSecondaryForMC 
            << std::noboolalpha << std::endl;
}

//____________________________________________________________________
//
// EOF
//
Commit	Line	Data
7984e5f7	1	//
	2	// This class calculates the exclusive charged particle density
	3	// in each for the 5 FMD rings.
	4	//
	5	// Input:
	6	// - 5 RingHistos objects - each with a number of vertex dependent
	7	// 2D histograms of the inclusive charge particle density
	8	//
	9	// Output:
	10	// - 5 RingHistos objects - each with a number of vertex dependent
	11	// 2D histograms of the exclusive charge particle density
	12	//
	13	// Corrections used:
	14	// - AliFMDCorrSecondaryMap;
	15	// - AliFMDCorrVertexBias
	16	// - AliFMDCorrMergingEfficiency
	17	//
72cc12cd	18	#include "AliFMDMCCorrector.h"
886dd429	19	#include <AliESDFMD.h>
	20	#include <TAxis.h>
	21	#include <TList.h>
	22	#include <TMath.h>
	23	#include "AliForwardCorrectionManager.h"
fb3430ac	24	#include "AliFMDCorrVertexBias.h"
886dd429	25	#include "AliLog.h"
	26	#include <TH2D.h>
	27	#include <TROOT.h>
	28	#include <TProfile2D.h>
e58000b7	29	#include <iostream>
72cc12cd	30	ClassImp(AliFMDMCCorrector)
886dd429	31	#if 0
	32	; // For Emacs
	33	#endif
	34
	35
	36	//____________________________________________________________________
72cc12cd	37	AliFMDMCCorrector::~AliFMDMCCorrector()
886dd429	38	{
7984e5f7	39	//
	40	// Destructor
	41	//
886dd429	42	if (fComps) fComps->Clear();
	43	if (fFMD1i) delete fFMD1i;
	44	if (fFMD2i) delete fFMD2i;
	45	if (fFMD2o) delete fFMD2o;
	46	if (fFMD3i) delete fFMD3i;
	47	if (fFMD3o) delete fFMD3o;
	48	}
	49
	50	//____________________________________________________________________
72cc12cd	51	AliFMDMCCorrector&
72cc12cd	52	AliFMDMCCorrector::operator=(const AliFMDMCCorrector& o)
886dd429	53	{
7984e5f7	54	//
	55	// Assignement operator
	56	//
	57	// Parameters:
	58	// o Object to assign from
	59	//
	60	// Return:
	61	// Reference to this object
	62	//
d015ecfe	63	if (&o == this) return *this;
72cc12cd	64	AliFMDCorrector::operator=(o);
5bb5d1f6	65	fSecondaryForMC = o.fSecondaryForMC;
886dd429	66	return *this;
	67	}
	68
	69	//____________________________________________________________________
	70	Bool_t
72cc12cd	71	AliFMDMCCorrector::CorrectMC(AliForwardUtil::Histos& hists,
886dd429	72	UShort_t vtxbin)
886dd429	73	{
7984e5f7	74	//
	75	// Do the calculations
	76	//
	77	// Parameters:
	78	// hists Cache of histograms
	79	// vtxBin Vertex bin
	80	//
	81	// Return:
	82	// true on successs
	83	//
5bb5d1f6	84	if ((!fUseSecondaryMap \|\| !fSecondaryForMC) && fUseVertexBias)
	85	return kTRUE;
	86
886dd429	87	AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
	88
	89	UShort_t uvb = vtxbin;
	90	for (UShort_t d=1; d<=3; d++) {
	91	UShort_t nr = (d == 1 ? 1 : 2);
	92	for (UShort_t q=0; q<nr; q++) {
	93	Char_t r = (q == 0 ? 'I' : 'O');
	94	TH2D* h = hists.Get(d,r);
7ec4d843	95
7ec4d843	96
5bb5d1f6	97	if (fUseSecondaryMap && fSecondaryForMC) {
7ec4d843	98	TH2D* bg = fcm.GetSecondaryMap()->GetCorrection(d,r,uvb);
	99	if (!bg) {
	100	AliWarning(Form("No secondary correction for FMDM%d%c in vertex bin %d",
	101	d, r, uvb));
	102	continue;
	103	}
	104	// Divide by primary/total ratio
	105	h->Divide(bg);
886dd429	106	}
7ec4d843	107	if (fUseVertexBias) {
7ec4d843	108	TH2D* ef = fcm.GetVertexBias()->GetCorrection(r, uvb);
e58000b7	109	if (!ef) {
7ec4d843	110	AliWarning(Form("No event vertex bias correction in vertex bin %d",
886dd429	111	uvb));
7ec4d843	112	continue;
	113	}
	114	// Divide by the event selection efficiency
	115	h->Divide(ef);
886dd429	116	}
886dd429	117	}
	118	}
	119
	120	return kTRUE;
	121	}
	122
	123	//____________________________________________________________________
	124	void
72cc12cd	125	AliFMDMCCorrector::Init(const TAxis& eAxis)
886dd429	126	{
7984e5f7	127	//
	128	// Initialize this object
	129	//
	130	// Parameters:
	131	// etaAxis Eta axis to use
	132	//
7ec4d843	133	AliFMDCorrector::Init(eAxis);
7ec4d843	134
886dd429	135	fFMD1i = Make(1,'I',eAxis);
	136	fFMD2i = Make(2,'I',eAxis);
	137	fFMD2o = Make(2,'O',eAxis);
	138	fFMD3i = Make(3,'I',eAxis);
	139	fFMD3o = Make(3,'O',eAxis);
	140
	141	fComps->Add(fFMD1i);
	142	fComps->Add(fFMD2i);
	143	fComps->Add(fFMD2o);
	144	fComps->Add(fFMD3i);
	145	fComps->Add(fFMD3o);
	146	}
	147
	148	//____________________________________________________________________
	149	TProfile2D*
72cc12cd	150	AliFMDMCCorrector::Make(UShort_t d, Char_t r,
886dd429	151	const TAxis& axis) const
886dd429	152	{
7984e5f7	153	//
	154	// MAke comparison profiles
	155	//
	156	// Parameters:
	157	// d Detector
	158	// r Ring
	159	// axis Eta axis
	160	//
	161	// Return:
	162	// Newly allocated profile object
	163	//
886dd429	164	TProfile2D* ret = new TProfile2D(Form("FMD%d%c_esd_vs_mc", d, r),
	165	Form("ESD/MC signal for FMD%d%c", d, r),
	166	axis.GetNbins(),
	167	axis.GetXmin(),
	168	axis.GetXmax(),
	169	(r == 'I' \|\| r == 'i') ? 20 : 40,
	170	0, 2*TMath::Pi());
	171	ret->GetXaxis()->SetTitle("#eta");
	172	ret->GetYaxis()->SetTitle("#varphi [degrees]");
	173	ret->GetZaxis()->SetTitle("#LT primary density ESD/MC#GT");
	174	ret->SetDirectory(0);
	175	return ret;
	176	}
	177	//____________________________________________________________________
	178	void
72cc12cd	179	AliFMDMCCorrector::Fill(UShort_t d, Char_t r, TH2* esd, TH2* mc)
886dd429	180	{
7984e5f7	181	//
	182	// Fill comparison profiles
	183	//
	184	// Parameters:
	185	// d Detector
	186	// r Ring
	187	// esd ESD histogram
	188	// mc MC histogram
	189	//
886dd429	190	if (!esd \|\| !mc) return;
	191	TProfile2D* p = 0;
	192	switch (d) {
	193	case 1: p = fFMD1i; break;
	194	case 2: p = (r == 'I' \|\| r == 'i' ? fFMD2i : fFMD2o); break;
	195	case 3: p = (r == 'I' \|\| r == 'i' ? fFMD3i : fFMD3o); break;
	196	}
	197	if (!p) return;
	198
	199	for (Int_t iEta = 1; iEta <= esd->GetNbinsX(); iEta++) {
	200	Double_t eta = esd->GetXaxis()->GetBinCenter(iEta);
	201	for (Int_t iPhi = 1; iPhi <= esd->GetNbinsY(); iPhi++) {
	202	Double_t phi = esd->GetYaxis()->GetBinCenter(iPhi);
	203	Double_t mEsd = esd->GetBinContent(iEta,iPhi);
	204	Double_t mMc = mc->GetBinContent(iEta,iPhi);
	205
	206	p->Fill(eta, phi, (mMc > 0 ? mEsd / mMc : 0));
	207	}
	208	}
	209	}
	210
	211	//____________________________________________________________________
	212	Bool_t
72cc12cd	213	AliFMDMCCorrector::CompareResults(AliForwardUtil::Histos& esd,
886dd429	214	AliForwardUtil::Histos& mc)
886dd429	215	{
7984e5f7	216	//
	217	// Compare the result of analysing the ESD for
	218	// the inclusive charged particle density to analysing
	219	// MC truth
	220	//
	221	// Parameters:
	222	// esd
	223	// mc
	224	//
	225	// Return:
	226	// true
	227	//
886dd429	228	Fill(1, 'I', esd.Get(1,'I'), mc.Get(1,'I'));
	229	Fill(2, 'I', esd.Get(2,'I'), mc.Get(2,'I'));
	230	Fill(2, 'O', esd.Get(2,'O'), mc.Get(2,'O'));
	231	Fill(3, 'I', esd.Get(3,'I'), mc.Get(3,'I'));
	232	Fill(3, 'O', esd.Get(3,'O'), mc.Get(3,'O'));
	233
	234	return kTRUE;
	235	}
	236
	237	//____________________________________________________________________
	238	void
72cc12cd	239	AliFMDMCCorrector::DefineOutput(TList* dir)
886dd429	240	{
7984e5f7	241	//
	242	// Output diagnostic histograms to directory
	243	//
	244	// Parameters:
	245	// dir List to write in
	246	//
72cc12cd	247	AliFMDCorrector::DefineOutput(dir);
886dd429	248	TList* d = static_cast<TList*>(dir->FindObject(GetName()));
	249
	250	fComps = new TList;
	251	fComps->SetName("esd_mc_comparison");
	252	d->Add(fComps);
	253	}
5bb5d1f6	254	//____________________________________________________________________
	255	void
	256	AliFMDMCCorrector::Print(Option_t* option) const
	257	{
	258	//
	259	// Print information
	260	// Parameters:
	261	// option Not used
	262	//
	263	char ind[gROOT->GetDirLevel()+1];
	264	for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
	265	ind[gROOT->GetDirLevel()] = '\0';
	266	AliFMDCorrector::Print(option);
	267	std::cout << std::boolalpha
	268	<< ind << " Use sec. map on MC: " << fSecondaryForMC
	269	<< std::noboolalpha << std::endl;
	270	}
886dd429	271
	272	//____________________________________________________________________
	273	//
	274	// EOF
	275	//
	276
	277
	278