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

#include "AliFMDMCHitEnergyFitter.h"
#include "AliMCEvent.h"
#include "AliESDEvent.h"
#include "AliStack.h"
#include "AliFMDHit.h"
#include "AliESDFMD.h"
#include "AliMCAuxHandler.h"
#include "AliGenEventHeader.h"
#include "AliLog.h"
#include "AliFMDEncodedEdx.h"
#include <TH1.h>
#include <TH2.h>
#include <TNtuple.h>
#include <TMath.h>


//____________________________________________________________________
AliFMDMCHitEnergyFitter::AliFMDMCHitEnergyFitter()
  : AliFMDEnergyFitter(),
    fSumPrimary(0),
    fSumSecondary(0),
    fIp(),
    fNTrack(0),
    fNPrimary(0),
    fTuple(0)
{
}

//____________________________________________________________________
AliFMDMCHitEnergyFitter::AliFMDMCHitEnergyFitter(const char* title,
						 Bool_t useTuple)
  : AliFMDEnergyFitter(title),
    fSumPrimary(0),
    fSumSecondary(0),
    fIp(),
    fNTrack(0),
    fNPrimary(0),
    fTuple(0)
{
  // Some defaults 
  fUseIncreasingBins = true;
  fDoFits            = true;
  fDoMakeObject      = false;
  fResidualMethod    = kNoResiduals;

  if (!useTuple) return;

  fTuple = new TNtuple("hits", "Hits", 
		       "primary:eta:phi:edep:length:"
		       "detector:ring:sector:strip:ipz");

}
//____________________________________________________________________
AliFMDMCHitEnergyFitter::~AliFMDMCHitEnergyFitter()
{}

//____________________________________________________________________
void
AliFMDMCHitEnergyFitter::CreateOutputObjects(TList* dir)
{
  DGUARD(fDebug, 2, "Create output objects w/MC hits (%p)", dir);
  AliFMDEnergyFitter::CreateOutputObjects(dir);
  // TList* d = static_cast<TList*>(l->FindObject(GetName()));
}
//____________________________________________________________________
Bool_t
AliFMDMCHitEnergyFitter::PreEvent(const AliMCEvent& mcInput)
{
  DGUARD(fDebug, 5, "Reset for event w/MC hits (%p)", &mcInput);
  fSumPrimary.Reset(0);
  fSumSecondary.Reset(0);

  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
AliFMDMCHitEnergyFitter::Event(const AliESDEvent& esdInput,
			       const AliMCEvent&  mcInput,
			       AliMCAuxHandler&   handler)
{
  DGUARD(fDebug, 5, "Process an event for MC hits (%p,%p,%p)",
	 &esdInput, &mcInput, &handler);
  PreEvent(mcInput);

  Int_t  nEntries = handler.GetNEntry();
  DMSG(fDebug,5, "We got a total of %d particles", nEntries);

  for (Int_t i = 0; i < nEntries; i++) {
    TClonesArray* hits = handler.GetEntryArray(i);
    if (!hits) continue;

    AccumulateHits(mcInput, *hits);
  }
  return PostEvent(esdInput);
}

//____________________________________________________________________
Bool_t
AliFMDMCHitEnergyFitter::AccumulateHits(const AliMCEvent& mcInput,
					const TClonesArray& hits)
{
  DGUARD(fDebug, 5, "Accumulate MC hit energy loss (%p,%p,%d)",
	 &mcInput, &hits, hits.GetEntries());

  AliStack*  stack = const_cast<AliMCEvent&>(mcInput).Stack();
  Float_t    cache[10];
  Int_t      nHit  = hits.GetEntries();
  Int_t      oldTr = -1;
  UShort_t   oldD  = 0;
  Char_t     oldR  = '\0';
  for (Int_t j = 0; j < nHit; j++) { 
	
    // Check the hit 
    AliFMDHit* hit  = static_cast<AliFMDHit*>(hits.At(j));
    if (!hit) continue;

    // Zero fill array 
    if (fTuple) for (Int_t k = 0; k < 10; k++) cache[k] = 0;

    // Get the track number 
    Int_t iTr = hit->Track();
    if (iTr < 0 || iTr >= fNTrack) 
      AliWarningF("Track # %d seems out of bounds [0,%d]", 
		  iTr, fNTrack-1);

    // Get the track 
    AliMCParticle* particle = 
      static_cast<AliMCParticle*>(mcInput.GetTrack(iTr));
	
    // Particle type 0: unknown, 1: primary, 2: secondary
    // Int_t flag = 0; - not used at the moment 

    // Check if this charged and a primary 
    if (particle) {
      // Only charged tracks 
      if (!particle->Charge() != 0) continue;
    }
    Bool_t   isPrimary = stack->IsPhysicalPrimary(iTr) && iTr < fNPrimary;
    Double_t eloss     = hit->Edep();
    Double_t length    = hit->Length();
    // Double_t eta  = particle->Eta();
    Double_t x         = hit->X() - fIp[0];
    Double_t y         = hit->Y() - fIp[1];
    Double_t z         = hit->Z() - fIp[2];
    Double_t r         = TMath::Sqrt(x*x + y*y);
    Double_t phi       = TMath::ATan2(y, x);
    Double_t theta     = TMath::ATan2(r, z);
    Double_t eta       = -TMath::Log(TMath::Tan(theta/2));
    UShort_t det       = hit->Detector();
    Char_t   rng       = hit->Ring();
    UShort_t sec       = hit->Sector();
    UShort_t str       = hit->Strip();

    Bool_t isDistinct = (iTr != oldTr) || (det != oldD) || (rng != oldR);
    if (isDistinct) {
      RingHistos* h = static_cast<RingHistos*>(GetRingHistos(det, rng));
      h->fKind->Fill(eta, isPrimary ? .5 : 1.5);
    }
    oldTr = iTr;
    oldD  = det;
    oldR  = rng;

    AliFMDFloatMap& sum     =  (isPrimary ? fSumPrimary : fSumSecondary);
    // Float_t         old     =  sum(det, rng, sec, str);
    sum(det, rng, sec, str) += eloss;
#if 0
    Info("", "FMD%d%c[%02d,%03d]-%s: Adding %10.7f to %10.7f -> %10.7f",
	 det, rng, sec, str, (isPrimary ? "1st" : "2nd"), 
	 eloss, old, sum(det, rng, sec, str));
#endif
    if (!fTuple) continue;

    // Leaves primary:eta:phi:edep:length:detector:ring:sector:strip:ipz
    cache[0] = isPrimary ? 1 : 0;
    cache[1] = eta;
    cache[2] = phi;
    cache[3] = eloss;
    cache[4] = length;
    cache[5] = hit->Detector();
    cache[6] = Int_t(hit->Ring());
    cache[7] = hit->Sector();
    cache[8] = hit->Strip();
    cache[9] = fIp[2];
    fTuple->Fill(cache);
  }
  return true;
}

//____________________________________________________________________
Bool_t AliFMDMCHitEnergyFitter::PostEvent(const AliESDEvent& input)
{
  DGUARD(fDebug, 5, "Convert MC hit energy loss (%p)", &input);

  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<AliFMDMCHitEnergyFitter::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 = fSumPrimary(d, r, s, t);
	  Float_t  totSec  = fSumSecondary(d, r, s, t);
	  Float_t  totAll  = totPrim + totSec;
	  Double_t esdEta  = esdFMD->Eta(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);
	}
      }
    }
  }
  return true;
}

//____________________________________________________________________
AliFMDEnergyFitter::RingHistos*
AliFMDMCHitEnergyFitter::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 AliFMDMCHitEnergyFitter::RingHistos(d,r);
}

//====================================================================
AliFMDMCHitEnergyFitter::RingHistos::RingHistos()
  : AliFMDEnergyFitter::RingHistos(), 
    fPrimary(0), 
    fSecondary(0),
    fKind(0)
{}

//____________________________________________________________________
AliFMDMCHitEnergyFitter::RingHistos::RingHistos(UShort_t d, Char_t r)
  : AliFMDEnergyFitter::RingHistos(d,r), 
    fPrimary(0),
    fSecondary(0),
    fKind(0)
{}
//____________________________________________________________________
TArrayD
AliFMDMCHitEnergyFitter::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
AliFMDMCHitEnergyFitter::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);
  fSecondary = Make("secondary","#sum#Delta_{true} of secondaries",
		    eAxis, maxDE, nDEbins, useIncrBin);
  fHist->SetXTitle("#Delta_{true}");
  fSecondary->SetXTitle("#Delta_{true}");
  fPrimary->SetXTitle("#Delta_{true}");

  fKind       = 0;
  if (eAxis.GetXbins()->GetArray()) 
    fKind = new TH2I("kind", "Particle types", 
		     eAxis.GetNbins(), eAxis.GetXbins()->GetArray(), 
		     2, 0, 2);
  else 
    fKind = new TH2I("kind", "Particle types", 
		     eAxis.GetNbins(), eAxis.GetXmin(), eAxis.GetXmax(),
		     2, 0, 2);
  fKind->SetXTitle("#eta");
  fKind->GetYaxis()->SetBinLabel(1, "Primary");
  fKind->GetYaxis()->SetBinLabel(2, "Secondary");
  fKind->SetDirectory(0);
  fKind->Sumw2();

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

//____________________________________________________________________
void
AliFMDMCHitEnergyFitter::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;
  }
}

//____________________________________________________________________
TObjArray*
AliFMDMCHitEnergyFitter::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("compartive");
  out->SetOwner();
  l->Add(out);
  return out;
}
Commit	Line	Data
0b7de667	1	#include "AliFMDMCHitEnergyFitter.h"
	2	#include "AliMCEvent.h"
	3	#include "AliESDEvent.h"
	4	#include "AliStack.h"
	5	#include "AliFMDHit.h"
	6	#include "AliESDFMD.h"
	7	#include "AliMCAuxHandler.h"
	8	#include "AliGenEventHeader.h"
	9	#include "AliLog.h"
	10	#include "AliFMDEncodedEdx.h"
	11	#include <TH1.h>
	12	#include <TH2.h>
	13	#include <TNtuple.h>
	14	#include <TMath.h>
	15
	16
	17	//____________________________________________________________________
	18	AliFMDMCHitEnergyFitter::AliFMDMCHitEnergyFitter()
	19	: AliFMDEnergyFitter(),
	20	fSumPrimary(0),
	21	fSumSecondary(0),
	22	fIp(),
	23	fNTrack(0),
	24	fNPrimary(0),
	25	fTuple(0)
	26	{
	27	}
	28
	29	//____________________________________________________________________
	30	AliFMDMCHitEnergyFitter::AliFMDMCHitEnergyFitter(const char* title,
	31	Bool_t useTuple)
	32	: AliFMDEnergyFitter(title),
	33	fSumPrimary(0),
	34	fSumSecondary(0),
	35	fIp(),
	36	fNTrack(0),
	37	fNPrimary(0),
	38	fTuple(0)
	39	{
	40	// Some defaults
	41	fUseIncreasingBins = true;
	42	fDoFits = true;
	43	fDoMakeObject = false;
	44	fResidualMethod = kNoResiduals;
	45
	46	if (!useTuple) return;
	47
	48	fTuple = new TNtuple("hits", "Hits",
	49	"primary:eta:phi:edep:length:"
	50	"detector:ring:sector:strip:ipz");
	51
	52	}
	53	//____________________________________________________________________
	54	AliFMDMCHitEnergyFitter::~AliFMDMCHitEnergyFitter()
	55	{}
	56
	57	//____________________________________________________________________
	58	void
	59	AliFMDMCHitEnergyFitter::CreateOutputObjects(TList* dir)
	60	{
	61	DGUARD(fDebug, 2, "Create output objects w/MC hits (%p)", dir);
	62	AliFMDEnergyFitter::CreateOutputObjects(dir);
	63	// TList* d = static_cast<TList*>(l->FindObject(GetName()));
	64	}
65	//____________________________________________________________________
66	Bool_t
67	AliFMDMCHitEnergyFitter::PreEvent(const AliMCEvent& mcInput)
68	{
69	DGUARD(fDebug, 5, "Reset for event w/MC hits (%p)", &mcInput);
70	fSumPrimary.Reset(0);
71	fSumSecondary.Reset(0);
72
73	AliMCEvent& mc = const_cast<AliMCEvent&>(mcInput);
74	AliHeader* header = mc.Header();
75	AliStack* stack = mc.Stack();
76	AliGenEventHeader* genHeader = header->GenEventHeader();
77
78	genHeader->PrimaryVertex(fIp);
79	fNTrack = stack->GetNtrack();
80	fNPrimary = stack->GetNprimary();
81
82	return true;
83	}
84
85	//____________________________________________________________________
86	Bool_t
87	AliFMDMCHitEnergyFitter::Event(const AliESDEvent& esdInput,
88	const AliMCEvent& mcInput,
89	AliMCAuxHandler& handler)
90	{
91	DGUARD(fDebug, 5, "Process an event for MC hits (%p,%p,%p)",
92	&esdInput, &mcInput, &handler);
93	PreEvent(mcInput);
94
95	Int_t nEntries = handler.GetNEntry();
96	DMSG(fDebug,5, "We got a total of %d particles", nEntries);
97
98	for (Int_t i = 0; i < nEntries; i++) {
99	TClonesArray* hits = handler.GetEntryArray(i);
100	if (!hits) continue;
101
102	AccumulateHits(mcInput, *hits);
103	}
104	return PostEvent(esdInput);
105	}
106
107	//____________________________________________________________________
108	Bool_t
109	AliFMDMCHitEnergyFitter::AccumulateHits(const AliMCEvent& mcInput,
110	const TClonesArray& hits)
111	{
112	DGUARD(fDebug, 5, "Accumulate MC hit energy loss (%p,%p,%d)",
113	&mcInput, &hits, hits.GetEntries());
114
115	AliStack* stack = const_cast<AliMCEvent&>(mcInput).Stack();
116	Float_t cache[10];
117	Int_t nHit = hits.GetEntries();
118	Int_t oldTr = -1;
119	UShort_t oldD = 0;
120	Char_t oldR = '\0';
121	for (Int_t j = 0; j < nHit; j++) {
122
123	// Check the hit
124	AliFMDHit* hit = static_cast<AliFMDHit*>(hits.At(j));
125	if (!hit) continue;
126
127	// Zero fill array
128	if (fTuple) for (Int_t k = 0; k < 10; k++) cache[k] = 0;
129
130	// Get the track number
131	Int_t iTr = hit->Track();
132	if (iTr < 0 \|\| iTr >= fNTrack)
133	AliWarningF("Track # %d seems out of bounds [0,%d]",
134	iTr, fNTrack-1);
135
136	// Get the track
137	AliMCParticle* particle =
138	static_cast<AliMCParticle*>(mcInput.GetTrack(iTr));
139
140	// Particle type 0: unknown, 1: primary, 2: secondary
141	// Int_t flag = 0; - not used at the moment
142
143	// Check if this charged and a primary
144	if (particle) {
145	// Only charged tracks
146	if (!particle->Charge() != 0) continue;
147	}
148	Bool_t isPrimary = stack->IsPhysicalPrimary(iTr) && iTr < fNPrimary;
149	Double_t eloss = hit->Edep();
150	Double_t length = hit->Length();
151	// Double_t eta = particle->Eta();
152	Double_t x = hit->X() - fIp[0];
153	Double_t y = hit->Y() - fIp[1];
154	Double_t z = hit->Z() - fIp[2];
155	Double_t r = TMath::Sqrt(xx + yy);
156	Double_t phi = TMath::ATan2(y, x);
157	Double_t theta = TMath::ATan2(r, z);
158	Double_t eta = -TMath::Log(TMath::Tan(theta/2));
159	UShort_t det = hit->Detector();
160	Char_t rng = hit->Ring();
161	UShort_t sec = hit->Sector();
162	UShort_t str = hit->Strip();
163
164	Bool_t isDistinct = (iTr != oldTr) \|\| (det != oldD) \|\| (rng != oldR);
165	if (isDistinct) {
166	RingHistos* h = static_cast<RingHistos*>(GetRingHistos(det, rng));
167	h->fKind->Fill(eta, isPrimary ? .5 : 1.5);
168	}
169	oldTr = iTr;
170	oldD = det;
171	oldR = rng;
172
173	AliFMDFloatMap& sum = (isPrimary ? fSumPrimary : fSumSecondary);
174	// Float_t old = sum(det, rng, sec, str);
175	sum(det, rng, sec, str) += eloss;
176	#if 0
177	Info("", "FMD%d%c[%02d,%03d]-%s: Adding %10.7f to %10.7f -> %10.7f",
178	det, rng, sec, str, (isPrimary ? "1st" : "2nd"),
179	eloss, old, sum(det, rng, sec, str));
180	#endif
181	if (!fTuple) continue;
182
183	// Leaves primary:eta:phi:edep:length:detector:ring:sector:strip:ipz
184	cache[0] = isPrimary ? 1 : 0;
185	cache[1] = eta;
186	cache[2] = phi;
187	cache[3] = eloss;
188	cache[4] = length;
189	cache[5] = hit->Detector();
190	cache[6] = Int_t(hit->Ring());
191	cache[7] = hit->Sector();
192	cache[8] = hit->Strip();
193	cache[9] = fIp[2];
194	fTuple->Fill(cache);
195	}
196	return true;
197	}
198
199	//____________________________________________________________________
200	Bool_t AliFMDMCHitEnergyFitter::PostEvent(const AliESDEvent& input)
201	{
202	DGUARD(fDebug, 5, "Convert MC hit energy loss (%p)", &input);
203
204	AliESDFMD* esdFMD = input.GetFMDData();
205	for (UShort_t d = 1; d <= 3; d++) {
206	UShort_t nQ = (d == 1 ? 1 : 2);
207	for (UShort_t q = 0; q < nQ; q++) {
208	UShort_t nS = (q == 0 ? 20 : 40);
209	UShort_t nT = (q == 0 ? 512 : 256);
210	Char_t r = (q == 0 ? 'I' : 'O');
211	RingHistos* h =
212	static_cast<AliFMDMCHitEnergyFitter::RingHistos*>(GetRingHistos(d,r));
213	if (!h) continue;
214
215	for (UShort_t s = 0; s < nS; s++) {
216	for (UShort_t t = 0; t < nT; t++) {
217	Float_t totPrim = fSumPrimary(d, r, s, t);
218	Float_t totSec = fSumSecondary(d, r, s, t);
219	Float_t totAll = totPrim + totSec;
220	Double_t esdEta = esdFMD->Eta(d, r, s, t);
221	if (totAll > 0) h->FillMC(0, esdEta, totAll);
222	if (totPrim > 0) h->FillMC(1, esdEta, totPrim);
223	if (totSec > 0) h->FillMC(2, esdEta, totSec);
224	}
225	}
226	}
227	}
228	return true;
229	}
230
231	//____________________________________________________________________
232	AliFMDEnergyFitter::RingHistos*
233	AliFMDMCHitEnergyFitter::CreateRingHistos(UShort_t d, Char_t r) const
234	{
235	// DGUARD(fDebug, 1, "Create Ring cache of MC hit energy loss (FMD%d%c)",d,r);
236	return new AliFMDMCHitEnergyFitter::RingHistos(d,r);
237	}
238
239	//====================================================================
240	AliFMDMCHitEnergyFitter::RingHistos::RingHistos()
241	: AliFMDEnergyFitter::RingHistos(),
242	fPrimary(0),
243	fSecondary(0),
244	fKind(0)
245	{}
246
247	//____________________________________________________________________
248	AliFMDMCHitEnergyFitter::RingHistos::RingHistos(UShort_t d, Char_t r)
249	: AliFMDEnergyFitter::RingHistos(d,r),
250	fPrimary(0),
251	fSecondary(0),
252	fKind(0)
253	{}
254	//____________________________________________________________________
255	TArrayD
256	AliFMDMCHitEnergyFitter::RingHistos::MakeIncreasingAxis(Int_t,
257	Double_t,
258	Double_t) const
259	{
260	// Make an increasing axis for ELoss distributions.
261	//
262	// We use the service function of AliFMDEncodedEdx to do this
263	TArrayD ret;
264	const AliFMDEncodedEdx::Spec& s = AliFMDEncodedEdx::GetdEAxis();
265	s.FillBinArray(ret);
266
267	return ret;
268	}
269	//____________________________________________________________________
270	void
271	AliFMDMCHitEnergyFitter::RingHistos::SetupForData(const TAxis& eAxis,
272	const TAxis& /cAxis/,
273	Double_t maxDE,
274	Int_t nDEbins,
275	Bool_t useIncrBin)
276	{
277	// AliFMDEnergyFitter::RingHistos::SetupForData(eAxis, cAxis, maxDE, nDEbins,
278	// useIncrBin);
279
280	fHist = Make("eloss", "#sum#Delta_{true} of all",
281	eAxis, maxDE, nDEbins, useIncrBin);
282	fPrimary = Make("primary", "#sum#Delta_{true} of primaries",
283	eAxis, maxDE, nDEbins, useIncrBin);
284	fSecondary = Make("secondary","#sum#Delta_{true} of secondaries",
285	eAxis, maxDE, nDEbins, useIncrBin);
286	fHist->SetXTitle("#Delta_{true}");
287	fSecondary->SetXTitle("#Delta_{true}");
288	fPrimary->SetXTitle("#Delta_{true}");
289
290	fKind = 0;
291	if (eAxis.GetXbins()->GetArray())
292	fKind = new TH2I("kind", "Particle types",
293	eAxis.GetNbins(), eAxis.GetXbins()->GetArray(),
294	2, 0, 2);
295	else
296	fKind = new TH2I("kind", "Particle types",
297	eAxis.GetNbins(), eAxis.GetXmin(), eAxis.GetXmax(),
298	2, 0, 2);
299	fKind->SetXTitle("#eta");
300	fKind->GetYaxis()->SetBinLabel(1, "Primary");
301	fKind->GetYaxis()->SetBinLabel(2, "Secondary");
302	fKind->SetDirectory(0);
303	fKind->Sumw2();
304
305	fList->Add(fHist);
306	fList->Add(fPrimary);
307	fList->Add(fSecondary);
308	fList->Add(fKind);
309	}
310
311	//____________________________________________________________________
312	void
313	AliFMDMCHitEnergyFitter::RingHistos::FillMC(UShort_t flag, Double_t eta,
314	Double_t mult)
315	{
316	switch (flag) {
317	// case 0: AliFMDEnergyFitter::RingHistos::Fill(false, eta, 0, mult); break;
318	case 0: fHist->Fill(eta, mult); break;
319	case 1: fPrimary->Fill(eta, mult); break;
320	case 2: fSecondary->Fill(eta, mult); break;
321	}
322	}
323
324	//____________________________________________________________________
325	TObjArray*
326	AliFMDMCHitEnergyFitter::RingHistos::Fit(TList* dir,
327	Double_t lowCut,
328	UShort_t nParticles,
329	UShort_t minEntries,
330	UShort_t minusBins,
331	Double_t relErrorCut,
332	Double_t chi2nuCut,
333	Double_t minWeight,
334	Double_t regCut,
335	EResidualMethod residuals) const
336	{
337	TObjArray* all = FitSlices(dir, "eloss", lowCut, nParticles, minEntries,
338	minusBins, relErrorCut, chi2nuCut, minWeight,
339	regCut, residuals, false, 0);
340	TObjArray* prim = FitSlices(dir, "primary", lowCut, nParticles, minEntries,
341	minusBins, relErrorCut, chi2nuCut, minWeight,
342	regCut, residuals, false, 0);
343	TObjArray* sec = FitSlices(dir, "secondary", lowCut, nParticles,
344	minEntries, minusBins, relErrorCut, chi2nuCut,
345	minWeight, regCut, residuals, false, 0);
346	if (!all \|\| !prim \|\| !sec) {
347	AliWarningF("Failed to get results for all(%p), primary(%p), and/or "
348	"secondary(%p)", all, prim, sec);
349	return 0;
350	}
351	// Already added to the sub-folders
352	// dir->Add(all);
353	// dir->Add(prim);
354	// dir->Add(sec);
355
356	Int_t nPrim = prim->GetEntriesFast();
357	TObjArray* out = new TObjArray;
358	for (Int_t i = 0; i < nPrim; i++) {
359	TH1* h = static_cast<TH1*>(prim->At(i));
360	if (!h) continue;
361
362	TAxis* xAxis = h->GetXaxis();
363	TH2* hh = 0;
364	if (xAxis->GetXbins()->GetArray())
365	hh = new TH2D(h->GetName(), h->GetTitle(),
366	xAxis->GetNbins(), xAxis->GetXbins()->GetArray(),
367	3, 0, 3);
368	else
369	hh = new TH2D(h->GetName(), h->GetTitle(),
370	xAxis->GetNbins(), xAxis->GetXmin(), xAxis->GetXmax(),
371	3, 0, 3);
372	hh->GetYaxis()->SetBinLabel(1, "Primary");
373	hh->GetYaxis()->SetBinLabel(2, "Secondary");
374	hh->GetYaxis()->SetBinLabel(3, "All");
375	hh->GetXaxis()->SetTitle("#eta");
376	out->Add(hh);
377	}
378	for (Int_t i = 0; i < nPrim; i++) {
379	TH2* h = static_cast<TH2*>(out->At(i));
380	if (!h) continue;
381
382	TH1* hp = static_cast<TH1*>(prim->At(i));
383	TH1* hs = static_cast<TH1*>(sec->At(i));
384	TH1* ha = static_cast<TH1*>(all->At(i));
385	TH1* hh[] = { hp, hs, ha };
386	for (Int_t j = 0; j < 3; j++) {
387	TH1* ph = hh[j];
388	if (!ph) continue;
389
390	for (Int_t k = 1; k <= ph->GetNbinsX(); k++) {
391	Double_t c = ph->GetBinContent(k);
392	Double_t e = ph->GetBinError(k);
393	h->SetBinContent(k, j+1, c);
394	h->SetBinError(k, j+1, e);
395	}
396	}
397	}
398	TList* l = GetOutputList(dir);
399	if (!l) return 0;
400
401	out->SetName("compartive");
402	out->SetOwner();
403	l->Add(out);
404	return out;
405	}