[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 
  //
  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, eta, 0, 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	//
	77	fFMD1i = Make(1,'I',eAxis);
	78	fFMD2i = Make(2,'I',eAxis);
	79	fFMD2o = Make(2,'O',eAxis);
	80	fFMD3i = Make(3,'I',eAxis);
	81	fFMD3o = Make(3,'O',eAxis);
	82	fFMD1iC = Make(1,'I');
	83	fFMD2iC = Make(2,'I');
	84	fFMD2oC = Make(2,'O');
	85	fFMD3iC = Make(3,'I');
	86	fFMD3oC = Make(3,'O');
	87
	88	fComps->Add(fFMD1i);
	89	fComps->Add(fFMD2i);
	90	fComps->Add(fFMD2o);
	91	fComps->Add(fFMD3i);
	92	fComps->Add(fFMD3o);
	93	fComps->Add(fFMD1iC);
	94	fComps->Add(fFMD2iC);
	95	fComps->Add(fFMD2oC);
	96	fComps->Add(fFMD3iC);
	97	fComps->Add(fFMD3oC);
	98	}
0bd4b00f	99
68d11054	100
0bd4b00f	101	//____________________________________________________________________
0bd4b00f	102	Bool_t
68d11054	103	AliFMDMCDensityCalculator::CalculateMC(const AliESDFMD& fmd,
68d11054	104	AliForwardUtil::Histos& hists)
0bd4b00f	105	{
7984e5f7	106	//
	107	// Calculate the charged particle density from the MC track references.
	108	//
	109	// Parameters:
	110	// event MC event
	111	// hists Histograms to fill
	112	// vz Interaction z coordinate @f$ v_z@f$
	113	// vtxBin bin corresponding to @f$ v_z@f$
	114	//
	115	// Return:
	116	// true on success
	117	//
68d11054	118	for (UShort_t d=1; d<=3; d++) {
	119	UShort_t nr = (d == 1 ? 1 : 2);
	120	for (UShort_t q=0; q<nr; q++) {
	121	Char_t r = (q == 0 ? 'I' : 'O');
	122	UShort_t ns= (q == 0 ? 20 : 40);
	123	UShort_t nt= (q == 0 ? 512 : 256);
	124	TH2D* h = hists.Get(d,r);
0bd4b00f	125
68d11054	126	for (UShort_t s=0; s<ns; s++) {
	127	for (UShort_t t=0; t<nt; t++) {
	128	Float_t mult = fmd.Multiplicity(d,r,s,t);
	129
	130	if (mult == 0 \|\| mult > 20) continue;
	131
	132	Float_t phi = fmd.Phi(d,r,s,t) / 180 * TMath::Pi();
	133	Float_t eta = fmd.Eta(d,r,s,t);
	134
	135	Float_t corr = Correction(d, r, s, t, eta, 0, false);
	136	Float_t sig = (corr <= 0 ? 0 : mult / corr);
	137	h->Fill(eta,phi,sig);
	138	}
	139	}
0bd4b00f	140	}
	141	}
	142	return kTRUE;
	143	}
	144
68d11054	145	//____________________________________________________________________
	146	TProfile2D*
	147	AliFMDMCDensityCalculator::Make(UShort_t d, Char_t r,
	148	const TAxis& axis) const
	149	{
7984e5f7	150	//
	151	// MAke comparison profiles
	152	//
	153	// Parameters:
	154	// d Detector
	155	// r Ring
	156	// axis Eta axis
	157	//
	158	// Return:
	159	// Newly allocated profile object
	160	//
68d11054	161	TProfile2D* ret = new TProfile2D(Form("FMD%d%c_esd_vs_mc", d, r),
	162	Form("ESD/MC signal for FMD%d%c", d, r),
	163	axis.GetNbins(),
	164	axis.GetXmin(),
	165	axis.GetXmax(),
	166	(r == 'I' \|\| r == 'i') ? 20 : 40,
	167	0, 2*TMath::Pi());
	168	ret->GetXaxis()->SetTitle("#eta");
	169	ret->GetYaxis()->SetTitle("#varphi [degrees]");
	170	ret->GetZaxis()->SetTitle("#LT incusive density ESD/MC#GT");
	171	ret->SetDirectory(0);
	172	return ret;
	173	}
	174	//____________________________________________________________________
	175	TH2D*
	176	AliFMDMCDensityCalculator::Make(UShort_t d, Char_t r) const
	177	{
7984e5f7	178	//
	179	// MAke comparison profiles
	180	//
	181	// Parameters:
	182	// d Detector
	183	// r Ring
	184	//
	185	// Return:
	186	// Newly allocated profile object
	187	//
68d11054	188	TH2D* ret = new TH2D(Form("FMD%d%c_corr_mc_esd", d, r),
	189	Form("ESD-MC correlation for FMD%d%c", d, r),
	190	200, 0, 20, 200, 0, 20);
	191	ret->GetXaxis()->SetTitle("m_{incl} (MC)");
	192	ret->GetYaxis()->SetTitle("#Delta/#Delta_{mp} (ESD)");
	193	ret->GetZaxis()->SetTitle("Correlation");
	194	ret->SetDirectory(0);
	195	return ret;
	196	}
	197	//____________________________________________________________________
	198	void
	199	AliFMDMCDensityCalculator::Fill(UShort_t d, Char_t r, TH2* esd, TH2* mc)
	200	{
7984e5f7	201	//
	202	// Fill comparison profiles
	203	//
	204	// Parameters:
	205	// d Detector
	206	// r Ring
	207	// esd ESD histogram
	208	// mc MC histogram
	209	//
68d11054	210	if (!esd \|\| !mc) return;
	211	TProfile2D* p = 0;
	212	TH2D* h = 0;
	213	switch (d) {
	214	case 1:
	215	p = fFMD1i;
	216	h = fFMD1iC;
	217	break;
	218	case 2:
	219	p = (r == 'I' \|\| r == 'i' ? fFMD2i : fFMD2o);
	220	h = (r == 'I' \|\| r == 'i' ? fFMD2iC : fFMD2oC);
	221	break;
	222	case 3:
	223	p = (r == 'I' \|\| r == 'i' ? fFMD3i : fFMD3o);
	224	h = (r == 'I' \|\| r == 'i' ? fFMD3iC : fFMD3oC);
	225	break;
	226	}
	227	if (!p) return;
	228
	229	for (Int_t iEta = 1; iEta <= esd->GetNbinsX(); iEta++) {
	230	Double_t eta = esd->GetXaxis()->GetBinCenter(iEta);
	231	for (Int_t iPhi = 1; iPhi <= esd->GetNbinsY(); iPhi++) {
	232	Double_t phi = esd->GetYaxis()->GetBinCenter(iPhi);
	233	Double_t mEsd = esd->GetBinContent(iEta,iPhi);
	234	Double_t mMc = mc->GetBinContent(iEta,iPhi);
	235
	236	p->Fill(eta, phi, (mMc > 0 ? mEsd / mMc : 0));
	237	h->Fill(mMc,mEsd);
	238	}
	239	}
	240	}
	241
0bd4b00f	242	//____________________________________________________________________
0bd4b00f	243	Bool_t
68d11054	244	AliFMDMCDensityCalculator::CompareResults(AliForwardUtil::Histos& esd,
68d11054	245	AliForwardUtil::Histos& mc)
0bd4b00f	246	{
7984e5f7	247	//
	248	// Compare the result of analysing the ESD for
	249	// the inclusive charged particle density to analysing
	250	// MC truth
	251	//
	252	// Parameters:
	253	// esd
	254	// mc
	255	//
	256	// Return:
	257	// true
	258	//
68d11054	259	Fill(1, 'I', esd.Get(1,'I'), mc.Get(1,'I'));
	260	Fill(2, 'I', esd.Get(2,'I'), mc.Get(2,'I'));
	261	Fill(2, 'O', esd.Get(2,'O'), mc.Get(2,'O'));
	262	Fill(3, 'I', esd.Get(3,'I'), mc.Get(3,'I'));
	263	Fill(3, 'O', esd.Get(3,'O'), mc.Get(3,'O'));
	264
	265	return kTRUE;
0bd4b00f	266	}
0bd4b00f	267
68d11054	268	//____________________________________________________________________
	269	void
	270	AliFMDMCDensityCalculator::DefineOutput(TList* dir)
	271	{
7984e5f7	272	//
	273	// Output diagnostic histograms to directory
	274	//
	275	// Parameters:
	276	// dir List to write in
	277	//
68d11054	278	AliFMDDensityCalculator::DefineOutput(dir);
	279	TList* d = static_cast<TList*>(dir->FindObject(GetName()));
	280
	281	fComps = new TList;
	282	fComps->SetName("esd_mc_comparison");
	283	d->Add(fComps);
	284
	285	}
0bd4b00f	286	//____________________________________________________________________
	287	//
	288	// EOF
	289	//
	290
	291
	292