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

// 
// This class calculates the inclusive charged particle density
// in each for the 5 FMD rings based on the MC truth.
//
// Input:
//   - AliMCEvent  MC truth event infromation
//
// Output:
//   - None
//
// Corrections used: 
//   - None
//
//
#include "AliFMDMCDensityCalculator.h"
#include <TMath.h>
#include "AliForwardCorrectionManager.h"
#include "AliFMDStripIndex.h"
#include "AliMCEvent.h"
#include "AliESDFMD.h"
#include "AliLog.h"
#include <TH2D.h>
#include <TProfile2D.h>

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

//____________________________________________________________________
AliFMDMCDensityCalculator::~AliFMDMCDensityCalculator()
{
  // 
  // 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;
  if (fFMD1iC) delete fFMD1iC;
  if (fFMD2iC) delete fFMD2iC;
  if (fFMD2oC) delete fFMD2oC;
  if (fFMD3iC) delete fFMD3iC;
  if (fFMD3oC) delete fFMD3oC;
}

//____________________________________________________________________
AliFMDMCDensityCalculator&
AliFMDMCDensityCalculator::operator=(const AliFMDMCDensityCalculator& o)
{
  // 
  // Assignement operator
  // 
  // Parameters:
  //    o Object to assign from 
  // 
  // Return:
  //    Reference to this object
  //
  AliFMDDensityCalculator::operator=(o);
  return *this;
}


//____________________________________________________________________
void
AliFMDMCDensityCalculator::Init(const TAxis& eAxis)
{
  // 
  // Initialize this object 
  // 
  // Parameters:
  //    etaAxis Eta axis to use 
  //
  AliFMDDensityCalculator::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); 
  fFMD1iC = Make(1,'I');
  fFMD2iC = Make(2,'I');
  fFMD2oC = Make(2,'O');
  fFMD3iC = Make(3,'I');
  fFMD3oC = Make(3,'O');
 
  fComps->Add(fFMD1i);
  fComps->Add(fFMD2i);
  fComps->Add(fFMD2o);
  fComps->Add(fFMD3i);
  fComps->Add(fFMD3o);
  fComps->Add(fFMD1iC);
  fComps->Add(fFMD2iC);
  fComps->Add(fFMD2oC);
  fComps->Add(fFMD3iC);
  fComps->Add(fFMD3oC);
}
    

//____________________________________________________________________
Bool_t
AliFMDMCDensityCalculator::CalculateMC(const AliESDFMD&        fmd,
				       AliForwardUtil::Histos& hists)
{
  // 
  // Calculate the charged particle density from the MC track references. 
  // 
  // Parameters:
  //    event  MC event
  //    hists  Histograms to fill
  //    vz     Interaction z coordinate @f$ v_z@f$
  //    vtxBin bin corresponding to @f$ v_z@f$
  // 
  // Return:
  //    true on success
  //
  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');
      UShort_t    ns= (q == 0 ?  20 :  40);
      UShort_t    nt= (q == 0 ? 512 : 256);
      TH2D*       h = hists.Get(d,r);

      for (UShort_t s=0; s<ns; s++) { 
	for (UShort_t t=0; t<nt; t++) {
	  Float_t mult = fmd.Multiplicity(d,r,s,t);
	  
	  if (mult == 0 || mult > 20) continue;

	  Float_t phi = fmd.Phi(d,r,s,t) / 180 * TMath::Pi();
	  Float_t eta = fmd.Eta(d,r,s,t);

	  Float_t corr = Correction(d, r, s, t, 0, eta, false);
	  Float_t sig  = (corr <= 0 ? 0 : mult / corr);
	  h->Fill(eta,phi,sig);
	}
      }
    }
  }
  return kTRUE;
}

//____________________________________________________________________
TProfile2D*
AliFMDMCDensityCalculator::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 incusive density ESD/MC#GT");
  ret->SetDirectory(0);
  return ret;
}
//____________________________________________________________________
TH2D*
AliFMDMCDensityCalculator::Make(UShort_t d, Char_t r) const
{
  // 
  // MAke comparison profiles
  // 
  // Parameters:
  //    d     Detector 
  //    r     Ring 
  // 
  // Return:
  //    Newly allocated profile object
  //
  TH2D* ret = new TH2D(Form("FMD%d%c_corr_mc_esd", d, r),
		       Form("ESD-MC correlation for FMD%d%c", d, r),
		       200, 0, 20, 200, 0, 20);
  ret->GetXaxis()->SetTitle("m_{incl} (MC)");
  ret->GetYaxis()->SetTitle("#Delta/#Delta_{mp} (ESD)");
  ret->GetZaxis()->SetTitle("Correlation");
  ret->SetDirectory(0);
  return ret;
}
//____________________________________________________________________
void
AliFMDMCDensityCalculator::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;
  TH2D*       h = 0;
  switch (d) { 
  case 1:  
    p = fFMD1i;                                   
    h = fFMD1iC;
    break;
  case 2:  
    p = (r == 'I' || r == 'i' ? fFMD2i  : fFMD2o); 
    h = (r == 'I' || r == 'i' ? fFMD2iC : fFMD2oC); 
    break;
  case 3:  
    p = (r == 'I' || r == 'i' ? fFMD3i  : fFMD3o); 
    h = (r == 'I' || r == 'i' ? fFMD3iC : fFMD3oC); 
    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));
      h->Fill(mMc,mEsd);
    }
  }
}

//____________________________________________________________________
Bool_t
AliFMDMCDensityCalculator::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
AliFMDMCDensityCalculator::DefineOutput(TList* dir)
{
  // 
  // Output diagnostic histograms to directory 
  // 
  // Parameters:
  //    dir List to write in
  //  
  AliFMDDensityCalculator::DefineOutput(dir);
  TList* d = static_cast<TList*>(dir->FindObject(GetName()));

  fComps = new TList;
  fComps->SetName("esd_mc_comparison");
  d->Add(fComps);

}
//____________________________________________________________________
//
// EOF
//
Commit	Line	Data
7984e5f7	1	//
	2	// This class calculates the inclusive charged particle density
	3	// in each for the 5 FMD rings based on the MC truth.
	4	//
	5	// Input:
	6	// - AliMCEvent MC truth event infromation
	7	//
	8	// Output:
	9	// - None
	10	//
	11	// Corrections used:
	12	// - None
	13	//
	14	//
0bd4b00f	15	#include "AliFMDMCDensityCalculator.h"
	16	#include <TMath.h>
	17	#include "AliForwardCorrectionManager.h"
	18	#include "AliFMDStripIndex.h"
	19	#include "AliMCEvent.h"
68d11054	20	#include "AliESDFMD.h"
0bd4b00f	21	#include "AliLog.h"
0bd4b00f	22	#include <TH2D.h>
68d11054	23	#include <TProfile2D.h>
0bd4b00f	24
	25	ClassImp(AliFMDMCDensityCalculator)
	26	#if 0
	27	; // For Emacs
	28	#endif
	29
68d11054	30	//____________________________________________________________________
	31	AliFMDMCDensityCalculator::~AliFMDMCDensityCalculator()
	32	{
7984e5f7	33	//
	34	// Destructor
	35	//
68d11054	36	if (fComps) fComps->Clear();
	37	if (fFMD1i) delete fFMD1i;
	38	if (fFMD2i) delete fFMD2i;
	39	if (fFMD2o) delete fFMD2o;
	40	if (fFMD3i) delete fFMD3i;
	41	if (fFMD3o) delete fFMD3o;
	42	if (fFMD1iC) delete fFMD1iC;
	43	if (fFMD2iC) delete fFMD2iC;
	44	if (fFMD2oC) delete fFMD2oC;
	45	if (fFMD3iC) delete fFMD3iC;
	46	if (fFMD3oC) delete fFMD3oC;
	47	}
0bd4b00f	48
	49	//____________________________________________________________________
	50	AliFMDMCDensityCalculator&
	51	AliFMDMCDensityCalculator::operator=(const AliFMDMCDensityCalculator& o)
	52	{
7984e5f7	53	//
	54	// Assignement operator
	55	//
	56	// Parameters:
	57	// o Object to assign from
	58	//
	59	// Return:
	60	// Reference to this object
	61	//
0bd4b00f	62	AliFMDDensityCalculator::operator=(o);
	63	return *this;
	64	}
	65
7984e5f7	66
	67	//____________________________________________________________________
	68	void
	69	AliFMDMCDensityCalculator::Init(const TAxis& eAxis)
	70	{
	71	//
	72	// Initialize this object
	73	//
	74	// Parameters:
	75	// etaAxis Eta axis to use
	76	//
7ec4d843	77	AliFMDDensityCalculator::Init(eAxis);
7984e5f7	78	fFMD1i = Make(1,'I',eAxis);
	79	fFMD2i = Make(2,'I',eAxis);
	80	fFMD2o = Make(2,'O',eAxis);
	81	fFMD3i = Make(3,'I',eAxis);
	82	fFMD3o = Make(3,'O',eAxis);
	83	fFMD1iC = Make(1,'I');
	84	fFMD2iC = Make(2,'I');
	85	fFMD2oC = Make(2,'O');
	86	fFMD3iC = Make(3,'I');
	87	fFMD3oC = Make(3,'O');
	88
	89	fComps->Add(fFMD1i);
	90	fComps->Add(fFMD2i);
	91	fComps->Add(fFMD2o);
	92	fComps->Add(fFMD3i);
	93	fComps->Add(fFMD3o);
	94	fComps->Add(fFMD1iC);
	95	fComps->Add(fFMD2iC);
	96	fComps->Add(fFMD2oC);
	97	fComps->Add(fFMD3iC);
	98	fComps->Add(fFMD3oC);
	99	}
0bd4b00f	100
68d11054	101
0bd4b00f	102	//____________________________________________________________________
0bd4b00f	103	Bool_t
68d11054	104	AliFMDMCDensityCalculator::CalculateMC(const AliESDFMD& fmd,
68d11054	105	AliForwardUtil::Histos& hists)
0bd4b00f	106	{
7984e5f7	107	//
	108	// Calculate the charged particle density from the MC track references.
	109	//
	110	// Parameters:
	111	// event MC event
	112	// hists Histograms to fill
	113	// vz Interaction z coordinate @f$ v_z@f$
	114	// vtxBin bin corresponding to @f$ v_z@f$
	115	//
	116	// Return:
	117	// true on success
	118	//
68d11054	119	for (UShort_t d=1; d<=3; d++) {
	120	UShort_t nr = (d == 1 ? 1 : 2);
	121	for (UShort_t q=0; q<nr; q++) {
	122	Char_t r = (q == 0 ? 'I' : 'O');
	123	UShort_t ns= (q == 0 ? 20 : 40);
	124	UShort_t nt= (q == 0 ? 512 : 256);
	125	TH2D* h = hists.Get(d,r);
0bd4b00f	126
68d11054	127	for (UShort_t s=0; s<ns; s++) {
	128	for (UShort_t t=0; t<nt; t++) {
	129	Float_t mult = fmd.Multiplicity(d,r,s,t);
	130
	131	if (mult == 0 \|\| mult > 20) continue;
	132
	133	Float_t phi = fmd.Phi(d,r,s,t) / 180 * TMath::Pi();
	134	Float_t eta = fmd.Eta(d,r,s,t);
	135
6feacd76	136	Float_t corr = Correction(d, r, s, t, 0, eta, false);
68d11054	137	Float_t sig = (corr <= 0 ? 0 : mult / corr);
	138	h->Fill(eta,phi,sig);
	139	}
	140	}
0bd4b00f	141	}
	142	}
	143	return kTRUE;
	144	}
	145
68d11054	146	//____________________________________________________________________
	147	TProfile2D*
	148	AliFMDMCDensityCalculator::Make(UShort_t d, Char_t r,
	149	const TAxis& axis) const
	150	{
7984e5f7	151	//
	152	// MAke comparison profiles
	153	//
	154	// Parameters:
	155	// d Detector
	156	// r Ring
	157	// axis Eta axis
	158	//
	159	// Return:
	160	// Newly allocated profile object
	161	//
68d11054	162	TProfile2D* ret = new TProfile2D(Form("FMD%d%c_esd_vs_mc", d, r),
	163	Form("ESD/MC signal for FMD%d%c", d, r),
	164	axis.GetNbins(),
	165	axis.GetXmin(),
	166	axis.GetXmax(),
	167	(r == 'I' \|\| r == 'i') ? 20 : 40,
	168	0, 2*TMath::Pi());
	169	ret->GetXaxis()->SetTitle("#eta");
	170	ret->GetYaxis()->SetTitle("#varphi [degrees]");
	171	ret->GetZaxis()->SetTitle("#LT incusive density ESD/MC#GT");
	172	ret->SetDirectory(0);
	173	return ret;
	174	}
	175	//____________________________________________________________________
	176	TH2D*
	177	AliFMDMCDensityCalculator::Make(UShort_t d, Char_t r) const
	178	{
7984e5f7	179	//
	180	// MAke comparison profiles
	181	//
	182	// Parameters:
	183	// d Detector
	184	// r Ring
	185	//
	186	// Return:
	187	// Newly allocated profile object
	188	//
68d11054	189	TH2D* ret = new TH2D(Form("FMD%d%c_corr_mc_esd", d, r),
	190	Form("ESD-MC correlation for FMD%d%c", d, r),
	191	200, 0, 20, 200, 0, 20);
	192	ret->GetXaxis()->SetTitle("m_{incl} (MC)");
	193	ret->GetYaxis()->SetTitle("#Delta/#Delta_{mp} (ESD)");
	194	ret->GetZaxis()->SetTitle("Correlation");
	195	ret->SetDirectory(0);
	196	return ret;
	197	}
	198	//____________________________________________________________________
	199	void
	200	AliFMDMCDensityCalculator::Fill(UShort_t d, Char_t r, TH2* esd, TH2* mc)
	201	{
7984e5f7	202	//
	203	// Fill comparison profiles
	204	//
	205	// Parameters:
	206	// d Detector
	207	// r Ring
	208	// esd ESD histogram
	209	// mc MC histogram
	210	//
68d11054	211	if (!esd \|\| !mc) return;
	212	TProfile2D* p = 0;
	213	TH2D* h = 0;
	214	switch (d) {
	215	case 1:
	216	p = fFMD1i;
	217	h = fFMD1iC;
	218	break;
	219	case 2:
	220	p = (r == 'I' \|\| r == 'i' ? fFMD2i : fFMD2o);
	221	h = (r == 'I' \|\| r == 'i' ? fFMD2iC : fFMD2oC);
	222	break;
	223	case 3:
	224	p = (r == 'I' \|\| r == 'i' ? fFMD3i : fFMD3o);
	225	h = (r == 'I' \|\| r == 'i' ? fFMD3iC : fFMD3oC);
	226	break;
	227	}
	228	if (!p) return;
	229
	230	for (Int_t iEta = 1; iEta <= esd->GetNbinsX(); iEta++) {
	231	Double_t eta = esd->GetXaxis()->GetBinCenter(iEta);
	232	for (Int_t iPhi = 1; iPhi <= esd->GetNbinsY(); iPhi++) {
	233	Double_t phi = esd->GetYaxis()->GetBinCenter(iPhi);
	234	Double_t mEsd = esd->GetBinContent(iEta,iPhi);
	235	Double_t mMc = mc->GetBinContent(iEta,iPhi);
	236
	237	p->Fill(eta, phi, (mMc > 0 ? mEsd / mMc : 0));
	238	h->Fill(mMc,mEsd);
	239	}
	240	}
	241	}
	242
0bd4b00f	243	//____________________________________________________________________
0bd4b00f	244	Bool_t
68d11054	245	AliFMDMCDensityCalculator::CompareResults(AliForwardUtil::Histos& esd,
68d11054	246	AliForwardUtil::Histos& mc)
0bd4b00f	247	{
7984e5f7	248	//
	249	// Compare the result of analysing the ESD for
	250	// the inclusive charged particle density to analysing
	251	// MC truth
	252	//
	253	// Parameters:
	254	// esd
	255	// mc
	256	//
	257	// Return:
	258	// true
	259	//
68d11054	260	Fill(1, 'I', esd.Get(1,'I'), mc.Get(1,'I'));
	261	Fill(2, 'I', esd.Get(2,'I'), mc.Get(2,'I'));
	262	Fill(2, 'O', esd.Get(2,'O'), mc.Get(2,'O'));
	263	Fill(3, 'I', esd.Get(3,'I'), mc.Get(3,'I'));
	264	Fill(3, 'O', esd.Get(3,'O'), mc.Get(3,'O'));
	265
	266	return kTRUE;
0bd4b00f	267	}
0bd4b00f	268
68d11054	269	//____________________________________________________________________
	270	void
	271	AliFMDMCDensityCalculator::DefineOutput(TList* dir)
	272	{
7984e5f7	273	//
	274	// Output diagnostic histograms to directory
	275	//
	276	// Parameters:
	277	// dir List to write in
	278	//
68d11054	279	AliFMDDensityCalculator::DefineOutput(dir);
	280	TList* d = static_cast<TList*>(dir->FindObject(GetName()));
	281
	282	fComps = new TList;
	283	fComps->SetName("esd_mc_comparison");
	284	d->Add(fComps);
	285
	286	}
0bd4b00f	287	//____________________________________________________________________
	288	//
	289	// EOF
	290	//
	291
	292
	293