[u/mrichter/AliRoot.git] / PWGLF / 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->Delete();
}

//____________________________________________________________________
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::SetupForData(const TAxis& eAxis)
{
  // 
  // Initialize this object 
  // 
  // Parameters:
  //    etaAxis Eta axis to use 
  //
  AliFMDDensityCalculator::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); 
  fFMD1iC = Make(1,'I');
  fFMD2iC = Make(2,'I');
  fFMD2oC = Make(2,'O');
  fFMD3iC = Make(3,'I');
  fFMD3oC = Make(3,'O');
  fFMD1iD = Make(1,'I', 20);
  fFMD2iD = Make(2,'I', 20);
  fFMD2oD = Make(2,'O', 20);
  fFMD3iD = Make(3,'I', 20);
  fFMD3oD = Make(3,'O', 20);
 
  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);
  fComps->Add(fFMD1iD);
  fComps->Add(fFMD2iD);
  fComps->Add(fFMD2oD);
  fComps->Add(fFMD3iD);
  fComps->Add(fFMD3oD);
}
    

//____________________________________________________________________
Bool_t
AliFMDMCDensityCalculator::CalculateMC(const AliESDFMD&        fmd,
				       AliForwardUtil::Histos& hists)
{
  // 
  // Calculate the charged particle density from the MC track references. 
  // 
  // Parameters:
  //    fmd    Forward 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, t, 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 [radians]");
  ret->GetZaxis()->SetTitle("#LT #deltaN_{ch,incl}=N_{ch,incl,ana}-N_{ch,incl,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),
		       41, -.5, 40.5, 410, -.5, 40.5);
  ret->GetXaxis()->SetTitle("m_{incl} (MC)");
  ret->GetYaxis()->SetTitle("m_{incl} from #Delta/#Delta_{mp} (ESD)");
  ret->GetZaxis()->SetTitle("Correlation");
  ret->SetDirectory(0);
  return ret;
}
//____________________________________________________________________
TH1D*
AliFMDMCDensityCalculator::Make(UShort_t d, Char_t r, Int_t max) const
{
  // 
  // MAke comparison profiles
  // 
  // Parameters:
  //    d     Detector 
  //    r     Ring 
  // 
  // Return:
  //    Newly allocated profile object
  //
  TH1D* ret = new TH1D(Form("FMD%d%c_diff_mc_esd", d, r),
		       Form("ESD-MC difference for FMD%d%c", d, r),
		       5*max,-max, max);
  ret->GetXaxis()->SetTitle("MC - ESD");
  ret->SetFillColor(AliForwardUtil::RingColor(d,r));
  ret->SetLineColor(AliForwardUtil::RingColor(d,r));
  ret->SetFillStyle(3001);
  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;
  TH1D*       x = 0;
  switch (d) { 
  case 1:  
    p = fFMD1i;                                   
    h = fFMD1iC;
    x = fFMD1iD;
    break;
  case 2:  
    p = (r == 'I' || r == 'i' ? fFMD2i  : fFMD2o); 
    h = (r == 'I' || r == 'i' ? fFMD2iC : fFMD2oC); 
    x = (r == 'I' || r == 'i' ? fFMD2iD : fFMD2oD); 
    break;
  case 3:  
    p = (r == 'I' || r == 'i' ? fFMD3i  : fFMD3o); 
    h = (r == 'I' || r == 'i' ? fFMD3iC : fFMD3oC); 
    x = (r == 'I' || r == 'i' ? fFMD3iD : fFMD3oD); 
    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);
      
      if (mMc <= 0 && mEsd <= 0) continue;
      x->Fill(mMc-mEsd);
      p->Fill(eta, phi, mEsd - mMc);
      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::CreateOutputObjects(TList* dir)
{
  // 
  // Output diagnostic histograms to directory 
  // 
  // Parameters:
  //    dir List to write in
  //  
  AliFMDDensityCalculator::CreateOutputObjects(dir);
  TList* d = static_cast<TList*>(dir->FindObject(GetName()));

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

}
//____________________________________________________________________
//
// EOF
//
Commit	Line	Data
	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	//
	15	#include "AliFMDMCDensityCalculator.h"
	16	#include <TMath.h>
	17	#include "AliForwardCorrectionManager.h"
	18	#include "AliFMDStripIndex.h"
	19	#include "AliMCEvent.h"
	20	#include "AliESDFMD.h"
	21	#include "AliLog.h"
	22	#include <TH2D.h>
	23	#include <TProfile2D.h>
	24
	25	ClassImp(AliFMDMCDensityCalculator)
	26	#if 0
	27	; // For Emacs
	28	#endif
	29
	30	//____________________________________________________________________
	31	AliFMDMCDensityCalculator::~AliFMDMCDensityCalculator()
	32	{
	33	//
	34	// Destructor
	35	//
	36	// if (fComps) fComps->Delete();
	37	}
	38
	39	//____________________________________________________________________
	40	AliFMDMCDensityCalculator&
	41	AliFMDMCDensityCalculator::operator=(const AliFMDMCDensityCalculator& o)
	42	{
	43	//
	44	// Assignement operator
	45	//
	46	// Parameters:
	47	// o Object to assign from
	48	//
	49	// Return:
	50	// Reference to this object
	51	//
	52	AliFMDDensityCalculator::operator=(o);
	53	return *this;
	54	}
	55
	56
	57	//____________________________________________________________________
	58	void
	59	AliFMDMCDensityCalculator::SetupForData(const TAxis& eAxis)
	60	{
	61	//
	62	// Initialize this object
	63	//
	64	// Parameters:
	65	// etaAxis Eta axis to use
	66	//
	67	AliFMDDensityCalculator::SetupForData(eAxis);
	68	fFMD1i = Make(1,'I',eAxis);
	69	fFMD2i = Make(2,'I',eAxis);
	70	fFMD2o = Make(2,'O',eAxis);
	71	fFMD3i = Make(3,'I',eAxis);
	72	fFMD3o = Make(3,'O',eAxis);
	73	fFMD1iC = Make(1,'I');
	74	fFMD2iC = Make(2,'I');
	75	fFMD2oC = Make(2,'O');
	76	fFMD3iC = Make(3,'I');
	77	fFMD3oC = Make(3,'O');
	78	fFMD1iD = Make(1,'I', 20);
	79	fFMD2iD = Make(2,'I', 20);
	80	fFMD2oD = Make(2,'O', 20);
	81	fFMD3iD = Make(3,'I', 20);
	82	fFMD3oD = Make(3,'O', 20);
	83
	84	fComps->Add(fFMD1i);
	85	fComps->Add(fFMD2i);
	86	fComps->Add(fFMD2o);
	87	fComps->Add(fFMD3i);
	88	fComps->Add(fFMD3o);
	89	fComps->Add(fFMD1iC);
	90	fComps->Add(fFMD2iC);
	91	fComps->Add(fFMD2oC);
	92	fComps->Add(fFMD3iC);
	93	fComps->Add(fFMD3oC);
	94	fComps->Add(fFMD1iD);
	95	fComps->Add(fFMD2iD);
	96	fComps->Add(fFMD2oD);
	97	fComps->Add(fFMD3iD);
	98	fComps->Add(fFMD3oD);
	99	}
	100
	101
	102	//____________________________________________________________________
	103	Bool_t
	104	AliFMDMCDensityCalculator::CalculateMC(const AliESDFMD& fmd,
	105	AliForwardUtil::Histos& hists)
	106	{
	107	//
	108	// Calculate the charged particle density from the MC track references.
	109	//
	110	// Parameters:
	111	// fmd Forward 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	//
	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);
	126
	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
	136	Float_t corr = Correction(d, r, t, eta, false);
	137	Float_t sig = (corr <= 0 ? 0 : mult / corr);
	138	h->Fill(eta,phi,sig);
	139	}
	140	}
	141	}
	142	}
	143	return kTRUE;
	144	}
	145
	146	//____________________________________________________________________
	147	TProfile2D*
	148	AliFMDMCDensityCalculator::Make(UShort_t d, Char_t r,
	149	const TAxis& axis) const
	150	{
	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	//
	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 [radians]");
	171	ret->GetZaxis()->SetTitle("#LT #deltaN_{ch,incl}=N_{ch,incl,ana}-N_{ch,incl,mc}#GT");
	172	ret->SetDirectory(0);
	173	return ret;
	174	}
	175	//____________________________________________________________________
	176	TH2D*
	177	AliFMDMCDensityCalculator::Make(UShort_t d, Char_t r) const
	178	{
	179	//
	180	// MAke comparison profiles
	181	//
	182	// Parameters:
	183	// d Detector
	184	// r Ring
	185	//
	186	// Return:
	187	// Newly allocated profile object
	188	//
	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	41, -.5, 40.5, 410, -.5, 40.5);
	192	ret->GetXaxis()->SetTitle("m_{incl} (MC)");
	193	ret->GetYaxis()->SetTitle("m_{incl} from #Delta/#Delta_{mp} (ESD)");
	194	ret->GetZaxis()->SetTitle("Correlation");
	195	ret->SetDirectory(0);
	196	return ret;
	197	}
	198	//____________________________________________________________________
	199	TH1D*
	200	AliFMDMCDensityCalculator::Make(UShort_t d, Char_t r, Int_t max) const
	201	{
	202	//
	203	// MAke comparison profiles
	204	//
	205	// Parameters:
	206	// d Detector
	207	// r Ring
	208	//
	209	// Return:
	210	// Newly allocated profile object
	211	//
	212	TH1D* ret = new TH1D(Form("FMD%d%c_diff_mc_esd", d, r),
	213	Form("ESD-MC difference for FMD%d%c", d, r),
	214	5*max,-max, max);
	215	ret->GetXaxis()->SetTitle("MC - ESD");
	216	ret->SetFillColor(AliForwardUtil::RingColor(d,r));
	217	ret->SetLineColor(AliForwardUtil::RingColor(d,r));
	218	ret->SetFillStyle(3001);
	219	ret->SetDirectory(0);
	220	return ret;
	221	}
	222	//____________________________________________________________________
	223	void
	224	AliFMDMCDensityCalculator::Fill(UShort_t d, Char_t r, TH2* esd, TH2* mc)
	225	{
	226	//
	227	// Fill comparison profiles
	228	//
	229	// Parameters:
	230	// d Detector
	231	// r Ring
	232	// esd ESD histogram
	233	// mc MC histogram
	234	//
	235	if (!esd \|\| !mc) return;
	236	TProfile2D* p = 0;
	237	TH2D* h = 0;
	238	TH1D* x = 0;
	239	switch (d) {
	240	case 1:
	241	p = fFMD1i;
	242	h = fFMD1iC;
	243	x = fFMD1iD;
	244	break;
	245	case 2:
	246	p = (r == 'I' \|\| r == 'i' ? fFMD2i : fFMD2o);
	247	h = (r == 'I' \|\| r == 'i' ? fFMD2iC : fFMD2oC);
	248	x = (r == 'I' \|\| r == 'i' ? fFMD2iD : fFMD2oD);
	249	break;
	250	case 3:
	251	p = (r == 'I' \|\| r == 'i' ? fFMD3i : fFMD3o);
	252	h = (r == 'I' \|\| r == 'i' ? fFMD3iC : fFMD3oC);
	253	x = (r == 'I' \|\| r == 'i' ? fFMD3iD : fFMD3oD);
	254	break;
	255	}
	256	if (!p) return;
	257
	258	for (Int_t iEta = 1; iEta <= esd->GetNbinsX(); iEta++) {
	259	Double_t eta = esd->GetXaxis()->GetBinCenter(iEta);
	260	for (Int_t iPhi = 1; iPhi <= esd->GetNbinsY(); iPhi++) {
	261	Double_t phi = esd->GetYaxis()->GetBinCenter(iPhi);
	262	Double_t mEsd = esd->GetBinContent(iEta,iPhi);
	263	Double_t mMc = mc->GetBinContent(iEta,iPhi);
	264
	265	if (mMc <= 0 && mEsd <= 0) continue;
	266	x->Fill(mMc-mEsd);
	267	p->Fill(eta, phi, mEsd - mMc);
	268	h->Fill(mMc,mEsd);
	269	}
	270	}
	271	}
	272
	273	//____________________________________________________________________
	274	Bool_t
	275	AliFMDMCDensityCalculator::CompareResults(AliForwardUtil::Histos& esd,
	276	AliForwardUtil::Histos& mc)
	277	{
	278	//
	279	// Compare the result of analysing the ESD for
	280	// the inclusive charged particle density to analysing
	281	// MC truth
	282	//
	283	// Parameters:
	284	// esd
	285	// mc
	286	//
	287	// Return:
	288	// true
	289	//
	290	Fill(1, 'I', esd.Get(1,'I'), mc.Get(1,'I'));
	291	Fill(2, 'I', esd.Get(2,'I'), mc.Get(2,'I'));
	292	Fill(2, 'O', esd.Get(2,'O'), mc.Get(2,'O'));
	293	Fill(3, 'I', esd.Get(3,'I'), mc.Get(3,'I'));
	294	Fill(3, 'O', esd.Get(3,'O'), mc.Get(3,'O'));
	295
	296	return kTRUE;
	297	}
	298
	299	//____________________________________________________________________
	300	void
	301	AliFMDMCDensityCalculator::CreateOutputObjects(TList* dir)
	302	{
	303	//
	304	// Output diagnostic histograms to directory
	305	//
	306	// Parameters:
	307	// dir List to write in
	308	//
	309	AliFMDDensityCalculator::CreateOutputObjects(dir);
	310	TList* d = static_cast<TList*>(dir->FindObject(GetName()));
	311
	312	fComps = new TList;
	313	fComps->SetOwner();
	314	fComps->SetName("esd_mc_comparison");
	315	d->Add(fComps);
	316
	317	}
	318	//____________________________________________________________________
	319	//
	320	// EOF
	321	//
	322
	323
	324