[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 "AliFMDCorrSecondaryMap.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::SetupForData(const TAxis& eAxis)
{
  // 
  // Initialize this object 
  // 
  // Parameters:
  //    etaAxis Eta axis to use 
  //
  AliFMDCorrector::SetupForData(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::CreateOutputObjects(TList* dir)
{
  // 
  // Output diagnostic histograms to directory 
  // 
  // Parameters:
  //    dir List to write in
  //  
  AliFMDCorrector::CreateOutputObjects(dir);
  TList* d = static_cast<TList*>(dir->FindObject(GetName()));

  fComps = new TList;
  fComps->SetName("esd_mc_comparison");
  d->Add(fComps);
}
#define PFB(N,FLAG)				\
  do {									\
    AliForwardUtil::PrintName(N);					\
    std::cout << std::boolalpha << (FLAG) << std::noboolalpha << std::endl; \
  } while(false)

//____________________________________________________________________
void
AliFMDMCCorrector::Print(Option_t* option) const
{
  // 
  // Print information
  // Parameters:
  //    option Not used 
  //
  AliFMDCorrector::Print(option);
  gROOT->IncreaseDirLevel();
  PFB("Use sec. map on MC", fSecondaryForMC);
  gROOT->DecreaseDirLevel();
}

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