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

#include "AliFMDMCTrackInspector.h"
#include "AliMCEvent.h"
#include "AliESDEvent.h"
#include "AliStack.h"
#include "AliESDFMD.h"
#include "AliGenEventHeader.h"
#include "AliHeader.h"
#include "AliLog.h"
#include "AliFMDEncodedEdx.h"
#include <TH1.h>
#include <TH2.h>
#include <TMath.h>
#include <TClonesArray.h>

//____________________________________________________________________
AliFMDMCTrackInspector::AliFMDMCTrackInspector()
  : AliFMDEnergyFitter(),
    fTracker(),
    fIp(),
    fNTrack(0),
    fNPrimary(0)
{
}

//____________________________________________________________________
AliFMDMCTrackInspector::AliFMDMCTrackInspector(const char* title)
  : AliFMDEnergyFitter(title),
    fTracker("tracker"),
    fIp(),
    fNTrack(0),
    fNPrimary(0)
{
  // Some defaults 
  fUseIncreasingBins = true;
  fDoFits            = true;
  fDoMakeObject      = false;
  fResidualMethod    = kNoResiduals;
}
//____________________________________________________________________
AliFMDMCTrackInspector::~AliFMDMCTrackInspector()
{}

//____________________________________________________________________
void
AliFMDMCTrackInspector::CreateOutputObjects(TList* dir)
{
  DGUARD(fDebug, 2, "Create output objects w/MC hits (%p)", dir);
  fTracker.CreateOutputObjects(dir);
  AliFMDEnergyFitter::CreateOutputObjects(dir);

  TIter    next(&fRingHistos);
  RingHistos* o = 0;
  while ((o = static_cast<RingHistos*>(next()))) {
    o->fBetaGammadEdx = static_cast<TH2*>(fTracker.GetBetaGammadEdx()->Clone());
    o->fBetaGammaEta  = static_cast<TH2*>(fTracker.GetBetaGammaEta()->Clone());
    o->fDedxEta       = static_cast<TH2*>(fTracker.GetDEdxEta()->Clone());
    o->fBetaGammadEdx->SetDirectory(0);
    o->fBetaGammaEta ->SetDirectory(0);
    o->fDedxEta      ->SetDirectory(0);
    o->fList->Add(o->fBetaGammadEdx);
    o->fList->Add(o->fBetaGammaEta );
    o->fList->Add(o->fDedxEta      );
  }

  // TList* d = static_cast<TList*>(l->FindObject(GetName()));
}
//____________________________________________________________________
Bool_t
AliFMDMCTrackInspector::PreEvent(const AliMCEvent& mcInput)
{
  DGUARD(fDebug, 5, "Reset for event w/MC hits (%p)", &mcInput);
  AliMCEvent&        mc        = const_cast<AliMCEvent&>(mcInput);
  AliHeader*         header    = mc.Header();
  AliStack*          stack     = mc.Stack();
  AliGenEventHeader* genHeader = header->GenEventHeader();
  
  genHeader->PrimaryVertex(fIp);
  fNTrack   = stack->GetNtrack();
  fNPrimary = stack->GetNprimary();

  return true;
}

//____________________________________________________________________
Bool_t
AliFMDMCTrackInspector::Event(const AliESDEvent& esdInput,
			      const AliMCEvent&  mcInput,
			      Double_t           cent)
{
  DGUARD(fDebug, 5, "Process an event for MC hits (%p,%p)",
	 &esdInput, &mcInput);
  PreEvent(mcInput);

  AliESDFMD*  esdFMD    = esdInput.GetFMDData();
  if (!esdFMD) return true;

  fTracker.Calculate(*esdFMD, mcInput, fIp[2], cent);

  return PostEvent();
}


//____________________________________________________________________
Bool_t AliFMDMCTrackInspector::PostEvent()
{
  DGUARD(fDebug, 5, "Fill MC hit energy loss");

  // AliESDFMD*  esdFMD    = input.GetFMDData();
  for (UShort_t d = 1; d <= 3; d++) { 
    UShort_t nQ = (d == 1 ? 1 : 2);
    for (UShort_t q = 0; q < nQ; q++) { 
      UShort_t    nS = (q == 0 ?  20 :  40);
      UShort_t    nT = (q == 0 ? 512 : 256);
      Char_t      r  = (q == 0 ? 'I' : 'O');
      RingHistos* h  = 
	static_cast<AliFMDMCTrackInspector::RingHistos*>(GetRingHistos(d,r));
      if (!h) continue;

      for (UShort_t s = 0; s < nS; s++) { 
	for (UShort_t t = 0; t < nT; t++) {
	  Float_t  totPrim = fTracker.GetPrimaries()(d, r, s, t);
	  Float_t  totSec  = fTracker.GetSecondaries()(d, r, s, t);
	  Float_t  totAll  = fTracker.GetAll()(d, r, s, t);
	  Double_t esdEta  = fTracker.GetEta()(d, r, s, t);
	  if (totAll  > 0) h->FillMC(0,  esdEta, totAll);
	  if (totPrim > 0) h->FillMC(1,  esdEta, totPrim);
	  if (totSec  > 0) h->FillMC(2,  esdEta, totSec);
	}
      }
    }
  }
  TClonesArray* hits = fTracker.GetHits();
  TIter next(hits);
  AliFMDMCTrackELoss::Hit* hit = 0;
  while ((hit = static_cast<AliFMDMCTrackELoss::Hit*>(next()))) {
    RingHistos* h  = 
      static_cast<AliFMDMCTrackInspector::RingHistos*>(GetRingHistos(hit->D(),
								     hit->R()));
    if (!h) continue;
    h->fBetaGammadEdx->Fill(hit->BetaGamma(), hit->DeDx());
    h->fBetaGammaEta->Fill(hit->Eta(), hit->BetaGamma());
    h->fDedxEta->Fill(hit->Eta(), hit->DeDx());
    
  }
  return true;
}

//____________________________________________________________________
AliFMDEnergyFitter::RingHistos*
AliFMDMCTrackInspector::CreateRingHistos(UShort_t d, Char_t r) const
{ 
  // DGUARD(fDebug, 1, "Create Ring cache of MC hit energy loss (FMD%d%c)",d,r);
  return new AliFMDMCTrackInspector::RingHistos(d,r);
}

//====================================================================
AliFMDMCTrackInspector::RingHistos::RingHistos()
  : AliFMDEnergyFitter::RingHistos(), 
    fPrimary(0), 
    fSecondary(0),
    fBetaGammadEdx(0),
    fBetaGammaEta(0),
    fDedxEta(0)
{}

//____________________________________________________________________
AliFMDMCTrackInspector::RingHistos::RingHistos(UShort_t d, Char_t r)
  : AliFMDEnergyFitter::RingHistos(d,r), 
    fPrimary(0),
    fSecondary(0),
    fBetaGammadEdx(0),
    fBetaGammaEta(0),
    fDedxEta(0)
{}
//____________________________________________________________________
TArrayD
AliFMDMCTrackInspector::RingHistos::MakeIncreasingAxis(Int_t, 
							Double_t,
							Double_t) const
{
  // Make an increasing axis for ELoss distributions. 
  // 
  // We use the service function of AliFMDEncodedEdx to do this 
  TArrayD ret;
  const AliFMDEncodedEdx::Spec& s = AliFMDEncodedEdx::GetdEAxis();
  s.FillBinArray(ret);

  return ret;
}
//____________________________________________________________________
void
AliFMDMCTrackInspector::RingHistos::SetupForData(const TAxis& eAxis, 
						  const TAxis& /*cAxis*/, 
						  Double_t maxDE, 
						  Int_t    nDEbins, 
						  Bool_t   useIncrBin)
{
  // AliFMDEnergyFitter::RingHistos::SetupForData(eAxis, cAxis, maxDE, nDEbins, 
  // useIncrBin);

  fHist      = Make("eloss", "#sum#Delta_{true} of all", 
		    eAxis, maxDE, nDEbins, useIncrBin);
  fPrimary   = Make("primary", "#sum#Delta_{true} of primaries", 
		    eAxis, maxDE, nDEbins, useIncrBin);
  fPrimary->SetMarkerStyle(24);
  fPrimary->SetMarkerSize(fPrimary->GetMarkerSize()*1.2);
  fSecondary = Make("secondary","#sum#Delta_{true} of secondaries",
 		    eAxis, maxDE, nDEbins, useIncrBin);
  fSecondary->SetMarkerStyle(25);
  fSecondary->SetMarkerSize(fSecondary->GetMarkerSize()*1.2);
  fHist->SetXTitle("#Delta_{true}");
  fSecondary->SetXTitle("#Delta_{true}");
  fPrimary->SetXTitle("#Delta_{true}");

  fList->Add(fHist);
  fList->Add(fPrimary);
  fList->Add(fSecondary);
}

//____________________________________________________________________
void
AliFMDMCTrackInspector::RingHistos::FillMC(UShort_t flag, Double_t eta, 
					    Double_t mult)
{
  switch (flag) { 
    // case 0: AliFMDEnergyFitter::RingHistos::Fill(false, eta, 0, mult); break;
  case 0: fHist->Fill(eta, mult); break;
  case 1: fPrimary->Fill(eta, mult); break;
  case 2: fSecondary->Fill(eta, mult); break;
  }
}

//____________________________________________________________________
void
AliFMDMCTrackInspector::RingHistos::Scale(TH1* dist) const
{
  // First scale to bin width 
  AliFMDEnergyFitter::RingHistos::Scale(dist);
  // Then smoothen the histogram 
  dist->Smooth(2);
}
  
//____________________________________________________________________
TObjArray*
AliFMDMCTrackInspector::RingHistos::Fit(TList*           dir, 
					Double_t         lowCut, 
					UShort_t         nParticles,
					UShort_t         minEntries,
					UShort_t         minusBins, 
					Double_t         relErrorCut, 
					Double_t         chi2nuCut,
					Double_t         minWeight,
					Double_t         regCut,
					EResidualMethod  residuals) const
{
  TObjArray* all  = FitSlices(dir, "eloss", lowCut, nParticles, minEntries, 
			      minusBins, relErrorCut, chi2nuCut, minWeight, 
			      regCut, residuals, false, 0);
  TObjArray* prim = FitSlices(dir, "primary", lowCut, nParticles, minEntries, 
			      minusBins, relErrorCut, chi2nuCut, minWeight, 
			      regCut, residuals, false, 0);
  TObjArray* sec  = FitSlices(dir, "secondary", lowCut, nParticles, 
			      minEntries, minusBins, relErrorCut, chi2nuCut, 
			      minWeight, regCut, residuals, false, 0);
  if (!all || !prim || !sec) {
    AliWarningF("Failed to get results for all(%p), primary(%p), and/or "
		"secondary(%p)", all, prim, sec);
    return 0;
  }
  // Already added to the sub-folders
  // dir->Add(all);
  // dir->Add(prim);
  // dir->Add(sec);

  Int_t      nPrim = prim->GetEntriesFast();
  TObjArray* out   = new TObjArray;
  for (Int_t i = 0; i < nPrim; i++) { 
    TH1* h = static_cast<TH1*>(prim->At(i));
    if (!h) continue;

    TAxis* xAxis  = h->GetXaxis();
    TH2*   hh     = 0;
    if (xAxis->GetXbins()->GetArray()) 
      hh = new TH2D(h->GetName(), h->GetTitle(), 
		    xAxis->GetNbins(), xAxis->GetXbins()->GetArray(),
		    3, 0, 3);
    else 
      hh = new TH2D(h->GetName(), h->GetTitle(), 
		    xAxis->GetNbins(), xAxis->GetXmin(), xAxis->GetXmax(),
		    3, 0, 3);
    hh->GetYaxis()->SetBinLabel(1, "Primary");
    hh->GetYaxis()->SetBinLabel(2, "Secondary");
    hh->GetYaxis()->SetBinLabel(3, "All");
    hh->GetXaxis()->SetTitle("#eta");
    out->Add(hh);
  }
  for (Int_t i = 0; i < nPrim; i++) { 
    TH2* h = static_cast<TH2*>(out->At(i));
    if (!h) continue;
    
    TH1* hp = static_cast<TH1*>(prim->At(i));
    TH1* hs = static_cast<TH1*>(sec->At(i));
    TH1* ha = static_cast<TH1*>(all->At(i));
    TH1* hh[] = { hp, hs, ha };
    for (Int_t j = 0; j < 3; j++) { 
      TH1* ph = hh[j];
      if (!ph) continue;

      for (Int_t k = 1; k <= ph->GetNbinsX(); k++) { 
	Double_t c = ph->GetBinContent(k);
	Double_t e = ph->GetBinError(k);
	h->SetBinContent(k, j+1, c);
	h->SetBinError(k, j+1, e);
      }
    }
  }
  TList* l = GetOutputList(dir);
  if (!l) return 0; 

  out->SetName("comparative");
  out->SetOwner();
  l->Add(out);
  return out;
}
// 
// EOF
//
Commit	Line	Data
7a3e34f4	1	#include "AliFMDMCTrackInspector.h"
	2	#include "AliMCEvent.h"
	3	#include "AliESDEvent.h"
	4	#include "AliStack.h"
	5	#include "AliESDFMD.h"
	6	#include "AliGenEventHeader.h"
	7	#include "AliHeader.h"
	8	#include "AliLog.h"
	9	#include "AliFMDEncodedEdx.h"
	10	#include <TH1.h>
	11	#include <TH2.h>
	12	#include <TMath.h>
	13	#include <TClonesArray.h>
	14
	15	//____________________________________________________________________
	16	AliFMDMCTrackInspector::AliFMDMCTrackInspector()
	17	: AliFMDEnergyFitter(),
	18	fTracker(),
	19	fIp(),
	20	fNTrack(0),
	21	fNPrimary(0)
	22	{
	23	}
	24
	25	//____________________________________________________________________
	26	AliFMDMCTrackInspector::AliFMDMCTrackInspector(const char* title)
	27	: AliFMDEnergyFitter(title),
	28	fTracker("tracker"),
	29	fIp(),
	30	fNTrack(0),
	31	fNPrimary(0)
	32	{
	33	// Some defaults
	34	fUseIncreasingBins = true;
	35	fDoFits = true;
	36	fDoMakeObject = false;
	37	fResidualMethod = kNoResiduals;
	38	}
	39	//____________________________________________________________________
	40	AliFMDMCTrackInspector::~AliFMDMCTrackInspector()
	41	{}
	42
	43	//____________________________________________________________________
	44	void
	45	AliFMDMCTrackInspector::CreateOutputObjects(TList* dir)
	46	{
	47	DGUARD(fDebug, 2, "Create output objects w/MC hits (%p)", dir);
	48	fTracker.CreateOutputObjects(dir);
	49	AliFMDEnergyFitter::CreateOutputObjects(dir);
	50
	51	TIter next(&fRingHistos);
	52	RingHistos* o = 0;
	53	while ((o = static_cast<RingHistos*>(next()))) {
	54	o->fBetaGammadEdx = static_cast<TH2*>(fTracker.GetBetaGammadEdx()->Clone());
	55	o->fBetaGammaEta = static_cast<TH2*>(fTracker.GetBetaGammaEta()->Clone());
	56	o->fDedxEta = static_cast<TH2*>(fTracker.GetDEdxEta()->Clone());
	57	o->fBetaGammadEdx->SetDirectory(0);
	58	o->fBetaGammaEta ->SetDirectory(0);
	59	o->fDedxEta ->SetDirectory(0);
	60	o->fList->Add(o->fBetaGammadEdx);
	61	o->fList->Add(o->fBetaGammaEta );
	62	o->fList->Add(o->fDedxEta );
	63	}
	64
65	// TList* d = static_cast<TList*>(l->FindObject(GetName()));
66	}
67	//____________________________________________________________________
68	Bool_t
69	AliFMDMCTrackInspector::PreEvent(const AliMCEvent& mcInput)
70	{
71	DGUARD(fDebug, 5, "Reset for event w/MC hits (%p)", &mcInput);
72	AliMCEvent& mc = const_cast<AliMCEvent&>(mcInput);
73	AliHeader* header = mc.Header();
74	AliStack* stack = mc.Stack();
75	AliGenEventHeader* genHeader = header->GenEventHeader();
76
77	genHeader->PrimaryVertex(fIp);
78	fNTrack = stack->GetNtrack();
79	fNPrimary = stack->GetNprimary();
80
81	return true;
82	}
83
84	//____________________________________________________________________
85	Bool_t
86	AliFMDMCTrackInspector::Event(const AliESDEvent& esdInput,
87	const AliMCEvent& mcInput,
88	Double_t cent)
89	{
90	DGUARD(fDebug, 5, "Process an event for MC hits (%p,%p)",
91	&esdInput, &mcInput);
92	PreEvent(mcInput);
93
94	AliESDFMD* esdFMD = esdInput.GetFMDData();
95	if (!esdFMD) return true;
96
97	fTracker.Calculate(*esdFMD, mcInput, fIp[2], cent);
98
99	return PostEvent();
100	}
101
102
103	//____________________________________________________________________
104	Bool_t AliFMDMCTrackInspector::PostEvent()
105	{
106	DGUARD(fDebug, 5, "Fill MC hit energy loss");
107
108	// AliESDFMD* esdFMD = input.GetFMDData();
109	for (UShort_t d = 1; d <= 3; d++) {
110	UShort_t nQ = (d == 1 ? 1 : 2);
111	for (UShort_t q = 0; q < nQ; q++) {
112	UShort_t nS = (q == 0 ? 20 : 40);
113	UShort_t nT = (q == 0 ? 512 : 256);
114	Char_t r = (q == 0 ? 'I' : 'O');
115	RingHistos* h =
116	static_cast<AliFMDMCTrackInspector::RingHistos*>(GetRingHistos(d,r));
117	if (!h) continue;
118
119	for (UShort_t s = 0; s < nS; s++) {
120	for (UShort_t t = 0; t < nT; t++) {
121	Float_t totPrim = fTracker.GetPrimaries()(d, r, s, t);
122	Float_t totSec = fTracker.GetSecondaries()(d, r, s, t);
123	Float_t totAll = fTracker.GetAll()(d, r, s, t);
124	Double_t esdEta = fTracker.GetEta()(d, r, s, t);
125	if (totAll > 0) h->FillMC(0, esdEta, totAll);
126	if (totPrim > 0) h->FillMC(1, esdEta, totPrim);
127	if (totSec > 0) h->FillMC(2, esdEta, totSec);
128	}
129	}
130	}
131	}
132	TClonesArray* hits = fTracker.GetHits();
133	TIter next(hits);
134	AliFMDMCTrackELoss::Hit* hit = 0;
135	while ((hit = static_cast<AliFMDMCTrackELoss::Hit*>(next()))) {
136	RingHistos* h =
137	static_cast<AliFMDMCTrackInspector::RingHistos*>(GetRingHistos(hit->D(),
138	hit->R()));
139	if (!h) continue;
140	h->fBetaGammadEdx->Fill(hit->BetaGamma(), hit->DeDx());
141	h->fBetaGammaEta->Fill(hit->Eta(), hit->BetaGamma());
142	h->fDedxEta->Fill(hit->Eta(), hit->DeDx());
143
144	}
145	return true;
146	}
147
148	//____________________________________________________________________
149	AliFMDEnergyFitter::RingHistos*
150	AliFMDMCTrackInspector::CreateRingHistos(UShort_t d, Char_t r) const
151	{
152	// DGUARD(fDebug, 1, "Create Ring cache of MC hit energy loss (FMD%d%c)",d,r);
153	return new AliFMDMCTrackInspector::RingHistos(d,r);
154	}
155
156	//====================================================================
157	AliFMDMCTrackInspector::RingHistos::RingHistos()
158	: AliFMDEnergyFitter::RingHistos(),
159	fPrimary(0),
160	fSecondary(0),
161	fBetaGammadEdx(0),
162	fBetaGammaEta(0),
163	fDedxEta(0)
164	{}
165
166	//____________________________________________________________________
167	AliFMDMCTrackInspector::RingHistos::RingHistos(UShort_t d, Char_t r)
168	: AliFMDEnergyFitter::RingHistos(d,r),
169	fPrimary(0),
170	fSecondary(0),
171	fBetaGammadEdx(0),
172	fBetaGammaEta(0),
173	fDedxEta(0)
174	{}
175	//____________________________________________________________________
176	TArrayD
177	AliFMDMCTrackInspector::RingHistos::MakeIncreasingAxis(Int_t,
178	Double_t,
179	Double_t) const
180	{
181	// Make an increasing axis for ELoss distributions.
182	//
183	// We use the service function of AliFMDEncodedEdx to do this
184	TArrayD ret;
185	const AliFMDEncodedEdx::Spec& s = AliFMDEncodedEdx::GetdEAxis();
186	s.FillBinArray(ret);
187
188	return ret;
189	}
190	//____________________________________________________________________
191	void
192	AliFMDMCTrackInspector::RingHistos::SetupForData(const TAxis& eAxis,
193	const TAxis& /cAxis/,
194	Double_t maxDE,
195	Int_t nDEbins,
196	Bool_t useIncrBin)
197	{
198	// AliFMDEnergyFitter::RingHistos::SetupForData(eAxis, cAxis, maxDE, nDEbins,
199	// useIncrBin);
200
201	fHist = Make("eloss", "#sum#Delta_{true} of all",
202	eAxis, maxDE, nDEbins, useIncrBin);
203	fPrimary = Make("primary", "#sum#Delta_{true} of primaries",
204	eAxis, maxDE, nDEbins, useIncrBin);
205	fPrimary->SetMarkerStyle(24);
206	fPrimary->SetMarkerSize(fPrimary->GetMarkerSize()*1.2);
207	fSecondary = Make("secondary","#sum#Delta_{true} of secondaries",
208	eAxis, maxDE, nDEbins, useIncrBin);
209	fSecondary->SetMarkerStyle(25);
210	fSecondary->SetMarkerSize(fSecondary->GetMarkerSize()*1.2);
211	fHist->SetXTitle("#Delta_{true}");
212	fSecondary->SetXTitle("#Delta_{true}");
213	fPrimary->SetXTitle("#Delta_{true}");
214
215	fList->Add(fHist);
216	fList->Add(fPrimary);
217	fList->Add(fSecondary);
218	}
219
220	//____________________________________________________________________
221	void
222	AliFMDMCTrackInspector::RingHistos::FillMC(UShort_t flag, Double_t eta,
223	Double_t mult)
224	{
225	switch (flag) {
226	// case 0: AliFMDEnergyFitter::RingHistos::Fill(false, eta, 0, mult); break;
227	case 0: fHist->Fill(eta, mult); break;
228	case 1: fPrimary->Fill(eta, mult); break;
229	case 2: fSecondary->Fill(eta, mult); break;
230	}
231	}
232
233	//____________________________________________________________________
92fdfa31	234	void
92fdfa31	235	AliFMDMCTrackInspector::RingHistos::Scale(TH1* dist) const
7a3e34f4	236	{
92fdfa31	237	// First scale to bin width
	238	AliFMDEnergyFitter::RingHistos::Scale(dist);
	239	// Then smoothen the histogram
	240	dist->Smooth(2);
7a3e34f4	241	}
92fdfa31	242
7a3e34f4	243	//____________________________________________________________________
	244	TObjArray*
	245	AliFMDMCTrackInspector::RingHistos::Fit(TList* dir,
	246	Double_t lowCut,
	247	UShort_t nParticles,
	248	UShort_t minEntries,
	249	UShort_t minusBins,
	250	Double_t relErrorCut,
	251	Double_t chi2nuCut,
	252	Double_t minWeight,
	253	Double_t regCut,
	254	EResidualMethod residuals) const
	255	{
7a3e34f4	256	TObjArray* all = FitSlices(dir, "eloss", lowCut, nParticles, minEntries,
	257	minusBins, relErrorCut, chi2nuCut, minWeight,
	258	regCut, residuals, false, 0);
	259	TObjArray* prim = FitSlices(dir, "primary", lowCut, nParticles, minEntries,
	260	minusBins, relErrorCut, chi2nuCut, minWeight,
	261	regCut, residuals, false, 0);
	262	TObjArray* sec = FitSlices(dir, "secondary", lowCut, nParticles,
	263	minEntries, minusBins, relErrorCut, chi2nuCut,
	264	minWeight, regCut, residuals, false, 0);
	265	if (!all \|\| !prim \|\| !sec) {
	266	AliWarningF("Failed to get results for all(%p), primary(%p), and/or "
	267	"secondary(%p)", all, prim, sec);
	268	return 0;
	269	}
	270	// Already added to the sub-folders
	271	// dir->Add(all);
	272	// dir->Add(prim);
	273	// dir->Add(sec);
	274
	275	Int_t nPrim = prim->GetEntriesFast();
	276	TObjArray* out = new TObjArray;
	277	for (Int_t i = 0; i < nPrim; i++) {
	278	TH1* h = static_cast<TH1*>(prim->At(i));
	279	if (!h) continue;
	280
	281	TAxis* xAxis = h->GetXaxis();
	282	TH2* hh = 0;
	283	if (xAxis->GetXbins()->GetArray())
	284	hh = new TH2D(h->GetName(), h->GetTitle(),
	285	xAxis->GetNbins(), xAxis->GetXbins()->GetArray(),
	286	3, 0, 3);
	287	else
	288	hh = new TH2D(h->GetName(), h->GetTitle(),
	289	xAxis->GetNbins(), xAxis->GetXmin(), xAxis->GetXmax(),
	290	3, 0, 3);
	291	hh->GetYaxis()->SetBinLabel(1, "Primary");
	292	hh->GetYaxis()->SetBinLabel(2, "Secondary");
	293	hh->GetYaxis()->SetBinLabel(3, "All");
	294	hh->GetXaxis()->SetTitle("#eta");
	295	out->Add(hh);
	296	}
	297	for (Int_t i = 0; i < nPrim; i++) {
	298	TH2* h = static_cast<TH2*>(out->At(i));
	299	if (!h) continue;
	300
	301	TH1* hp = static_cast<TH1*>(prim->At(i));
	302	TH1* hs = static_cast<TH1*>(sec->At(i));
	303	TH1* ha = static_cast<TH1*>(all->At(i));
	304	TH1* hh[] = { hp, hs, ha };
	305	for (Int_t j = 0; j < 3; j++) {
	306	TH1* ph = hh[j];
	307	if (!ph) continue;
	308
	309	for (Int_t k = 1; k <= ph->GetNbinsX(); k++) {
	310	Double_t c = ph->GetBinContent(k);
	311	Double_t e = ph->GetBinError(k);
	312	h->SetBinContent(k, j+1, c);
	313	h->SetBinError(k, j+1, e);
	314	}
	315	}
	316	}
	317	TList* l = GetOutputList(dir);
	318	if (!l) return 0;
	319
320	out->SetName("comparative");
321	out->SetOwner();
322	l->Add(out);
323	return out;
324	}
325	//
326	// EOF
327	//