[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TOF / PbPb276 / macros / HistoUtils.C

TH1 *
HistoUtils_weightedmean(TH1 *h1, TH1 *h2)
{

  TH1 *ho = (TH1 *)h1->Clone("ho");
  ho->Reset();
  Double_t val1, val2, w1, w2, mean, meane;
  for (Int_t i = 0; i < ho->GetNbinsX(); i++) {
    val1 = h1->GetBinContent(i + 1);
    w1 = 1. / (h1->GetBinError(i + 1) * h1->GetBinError(i + 1));
    val2 = h2->GetBinContent(i + 1);
    w2 = 1. / (h2->GetBinError(i + 1) * h2->GetBinError(i + 1));

    if (val1 == 0 && val2 == 0) continue;

    mean = (w1 * val1 + w2 * val2) / (w1 + w2);
    meane = TMath::Sqrt(1. / (w1 + w2));

    ho->SetBinContent(i + 1, mean);
    ho->SetBinError(i + 1, meane);
  }

  return ho;
}

//__________________________________________________________________

TH1 *
HistoUtils_smartdifference(TH1 *hnum, TH1 *hden)
{

  TH1 *hr = (TH1 *)hnum->Clone("hr");
  hr->Reset();
  Double_t ref;
  for (Int_t i = 0; i < hr->GetNbinsX(); i++) {
    if (hnum->GetBinError(i + 1) <= 0.) continue;
    ref = hden->Interpolate(hr->GetBinCenter(i + 1));
    if (ref <= 0.) continue;
    hr->SetBinContent(i + 1, (hnum->GetBinContent(i + 1) - ref) / ref);
    hr->SetBinError(i + 1, hnum->GetBinError(i + 1) / ref);
  }
  return hr;
}

//__________________________________________________________________

TH1 *
HistoUtils_smartratio(TH1 *hnum, TH1 *hden)
{

  TH1 *hr = (TH1 *)hnum->Clone("hr");
  hr->Reset();
  Double_t ref;
  for (Int_t i = 0; i < hr->GetNbinsX(); i++) {
    if (hnum->GetBinError(i + 1) <= 0.) continue;
    ref = hden->Interpolate(hr->GetBinCenter(i + 1));
    if (ref <= 0.) continue;
    hr->SetBinContent(i + 1, hnum->GetBinContent(i + 1) / ref);
    hr->SetBinError(i + 1, hnum->GetBinError(i + 1) / ref);
  }
  return hr;
}

//__________________________________________________________________

TH1 *
HistoUtils_smartratio(TH1 *hnum, TGraph *hden)
{

  TH1 *hr = (TH1 *)hnum->Clone("hr");
  hr->Reset();
  Double_t ref;
  for (Int_t i = 0; i < hr->GetNbinsX(); i++) {
    if (hnum->GetBinError(i + 1) <= 0.) continue;
    ref = hden->Eval(hr->GetBinCenter(i + 1));
    if (ref <= 0.) continue;
    hr->SetBinContent(i + 1, hnum->GetBinContent(i + 1) / ref);
    hr->SetBinError(i + 1, hnum->GetBinError(i + 1) / ref);
  }
  return hr;
}

//__________________________________________________________________

HistoUtils_drawthemall(const Char_t *filename, Int_t sleepms = 100)
{
  TFile *f1 = TFile::Open(filename);
  TList *l1 = f1->GetListOfKeys();
  TObject *o;
  Char_t *name;
  for (Int_t i = 0; i < l1->GetEntries(); i++) {
    name = l1->At(i)->GetName();
    o = f1->Get(name);
    o->Draw();
    gPad->Update();
    gSystem->Sleep(sleepms);
  }
  f1->Close();
}

//__________________________________________________________________

HistoUtils_autoratio(const Char_t *f1name, const Char_t *f2name, const Char_t *outname, const Char_t *title = NULL)
{

  TFile *f1 = TFile::Open(f1name);
  TFile *f2 = TFile::Open(f2name);
  TFile *fo = TFile::Open(outname, "RECREATE");
  TList *l1 = f1->GetListOfKeys();


  TH1D *h1, *h2, *hr;
  Char_t *name;
  for (Int_t i = 0; i < l1->GetEntries(); i++) {
    name = l1->At(i)->GetName();
    h1 = (TH1D *)f1->Get(name);
    h2 = (TH1D *)f2->Get(name);
    if (!h1 || !h2) continue;
    hr = new TH1D(*h1);
    hr->Divide(h2);
    if (title)
      hr->SetTitle(title);
    fo->cd();
    hr->Write();
  }

  delete hr;
  f1->Close();
  f2->Close();
  fo->Close();

}

//__________________________________________________________________

HistoUtils_autosystematics(const Char_t *f1name, const Char_t *f2name, const Char_t *outname, const Char_t *title = NULL)
{

  TFile *f1 = TFile::Open(f1name);
  TFile *f2 = TFile::Open(f2name);
  if (!f1 || !f1->IsOpen() || !f2 || !f2->IsOpen()) return;
  TFile *fo = TFile::Open(outname, "RECREATE");
  TList *l1 = f1->GetListOfKeys();


  TH1D *h1, *h2, *hd;
  Char_t *name;
  for (Int_t i = 0; i < l1->GetEntries(); i++) {
    name = l1->At(i)->GetName();
    h1 = (TH1D *)f1->Get(name);
    h2 = (TH1D *)f2->Get(name);
    if (!h1 || !h2) continue;

    hd = new TH1D(*h1);
    Double_t val1, val2, vald, vale1, vale2, valde;
    for (Int_t ii = 0; ii < hd->GetNbinsX(); ii++) {
      val1 = h1->GetBinContent(ii + 1);
      vale1 = h1->GetBinError(ii + 1);
      val2 = h2->GetBinContent(ii + 1);
      vale2 = h2->GetBinError(ii + 1);
      if (val2 == 0.) continue;
      vald = (val1 - val2) / val2;
      valde = TMath::Sqrt(TMath::Abs((vale1 * vale1 - vale2 * vale2))) / val2;
      hd->SetBinContent(ii + 1, vald);
      hd->SetBinError(ii + 1, valde);
    }

    if (title)
      hd->SetTitle(title);
    fo->cd();
    hd->Write();
    delete hd;
  }

  f1->Close();
  f2->Close();
  fo->Close();

}

//__________________________________________________________________

TH1 *
HistoUtils_ratio(const Char_t *f1name, const Char_t *f2name, const Char_t *h1name, const Char_t *h2name)
{

  if (!f2name) f2name = f1name;
  TFile *f1 = TFile::Open(f1name);
  TFile *f2 = TFile::Open(f2name);

  if (!h2name) h2name = h1name;
  TH1 *h1 = (TH1 *)f1->Get(h1name);
  TH1 *h2 = (TH1 *)f2->Get(h2name);

  TH1 *hr = h1->Clone("hr");
  hr->Sumw2();
  hr->Divide(h2);

  return hr;
  
}

//__________________________________________________________________

TH1 *
HistoUtils_systematics(const Char_t *f1name, const Char_t *f2name, const Char_t *h1name, const Char_t *h2name)
{

  if (!f2name) f2name = f1name;
  TFile *f1 = TFile::Open(f1name);
  TFile *f2 = TFile::Open(f2name);

  if (!h2name) h2name = h1name;
  TH1 *h1 = (TH1 *)f1->Get(h1name);
  TH1 *h2 = (TH1 *)f2->Get(h2name);

  TH1 *hd = h1->Clone("hd");
  Double_t val1, val2, vald, vale1, vale2, valde;
  for (Int_t i = 0; i < h1->GetNbinsX(); i++) {
    val1 = h1->GetBinContent(i + 1);
    vale1 = h1->GetBinError(i + 1);
    val2 = h2->GetBinContent(i + 1);
    vale2 = h2->GetBinError(i + 1);
    if (val2 == 0.) continue;
    vald = (val1 - val2) / val2;
    valde = TMath::Sqrt(TMath::Abs((vale1 * vale1 - vale2 * vale2))) / val2;
    hd->SetBinContent(i + 1, vald);
    hd->SetBinError(i + 1, valde);
  }

  return hd;
  
}

//__________________________________________________________________

HistoUtils_BinLogX(TH1 *h)
{
  TAxis *axis = h->GetXaxis();
  Int_t bins = axis->GetNbins();
  Axis_t from = axis->GetXmin();
  Axis_t to = axis->GetXmax();
  Axis_t width = (to - from) / bins;
  Axis_t *new_bins = new Axis_t[bins + 1];
  for (int i = 0; i <= bins; i++) new_bins[i] = TMath::Power(10, from + i * width);
  axis->Set(bins, new_bins);
  delete new_bins;
}

//__________________________________________________________________

HistoUtils_BinNormX(TH1 *h)
{
  TAxis *axis = h->GetXaxis();
  Int_t bins = axis->GetNbins();
  Double_t c, ec;
  Double_t w;
  for (Int_t i = 0; i < bins; i++) {
    c = h->GetBinContent(i + 1);
    ec = h->GetBinError(i + 1);
    w = axis->GetBinWidth(i + 1);
    c /= w;
    ec /= w;
    h->SetBinContent(i + 1, c);
    h->SetBinError(i + 1, ec);
  }
}

//__________________________________________________________________

TObjArray *
HistoUtils_FitPeak(TF1 *fitFunc, TH2 *h, Float_t startSigma, Float_t nSigmaMin, Float_t nSigmaMax, Int_t minIntegral = 100., const Char_t *basename = "hParam", Bool_t monitor = kFALSE)
{

  /* gaus function if not specified */
  if (!fitFunc)
    fitFunc = (TF1*)gROOT->GetFunction("gaus");

  /* prepare output histos */
  TObjArray *outArray = new TObjArray();
  Int_t npars = fitFunc->GetNpar();
  TH1D *hpx = h->ProjectionX("hpx");
  TH1D *hParam;
  for (Int_t ipar = 0; ipar < npars; ipar++) {
    hParam = new TH1D(*hpx);
    hParam->SetName(Form("%s_%d", basename, ipar));
    hParam->Reset();
    outArray->Add(hParam);
  }

  /* loop over x-bins */
  for (Int_t ibin = 0; ibin < hpx->GetNbinsX(); ibin++) {

    /* check integral */
    if (hpx->GetBinContent(ibin + 1) < minIntegral) continue;
    /* projection y */
    TH1D *hpy = h->ProjectionY("hpy", ibin + 1, ibin + 1);
    /* fit peak */
    if (HistoUtils_FitPeak(fitFunc, hpy, startSigma, nSigmaMin, nSigmaMax) != 0) {
      delete hpy;
      continue;
    }
    /* setup output histos */
    for (Int_t ipar = 0; ipar < npars; ipar++) {
      hParam = (TH1D *)outArray->At(ipar);
      hParam->SetBinContent(ibin + 1, fitFunc->GetParameter(ipar));
      hParam->SetBinError(ibin + 1, fitFunc->GetParError(ipar));
    }
    /* monitor */
    if (monitor) {
      hpy->SetMarkerStyle(20);
      hpy->SetMarkerColor(4);
      hpy->Draw("E1");
      fitFunc->Draw("same");
      gPad->Update();
      getchar();
    }
    /* delete */
    delete hpy;

  }
  /* delete */
  delete hpx;
  /* return output array */
  return outArray;
}

//__________________________________________________________________

Int_t
HistoUtils_FitPeak(TF1 *fitFunc, TH1 *h, Float_t startSigma, Float_t nSigmaMin, Float_t nSigmaMax)
{

  Double_t fitCent = h->GetBinCenter(h->GetMaximumBin());
  Double_t fitMin = fitCent - nSigmaMin * startSigma;
  Double_t fitMax = fitCent + nSigmaMax * startSigma;
  if (fitMin < h->GetXaxis()->GetXmin()) fitMin = h->GetXaxis()->GetXmin();
  if (fitMax > h->GetXaxis()->GetXmax()) fitMax = h->GetXaxis()->GetXmax();
  fitFunc->SetParameter(1, fitCent);
  fitFunc->SetParameter(2, startSigma);
  Int_t fitres = h->Fit(fitFunc, "WWq0", "", fitMin, fitMax);
  if (fitres != 0) return fitres;
  /* refit with better range */
  for (Int_t i = 0; i < 3; i++) {
    fitCent = fitFunc->GetParameter(1);
    fitMin = fitCent - nSigmaMin * fitFunc->GetParameter(2);
    fitMax = fitCent + nSigmaMax * fitFunc->GetParameter(2);
    if (fitMin < h->GetXaxis()->GetXmin()) fitMin = h->GetXaxis()->GetXmin();
    if (fitMax > h->GetXaxis()->GetXmax()) fitMax = h->GetXaxis()->GetXmax();
    fitres = h->Fit(fitFunc, "q0", "", fitMin, fitMax);
    if (fitres != 0) return fitres;
  }
  return fitres;

}

//__________________________________________________________________

void
HistoUtils_Function2Profile(TF1 *fin, TProfile *p, Int_t ntry = 10000)
{

  TF1 *f = new TF1(*fin);
  Int_t npars = f->GetNpar();
  Double_t par[1000];
  Double_t pare[1000];
  for (Int_t ipar = 0; ipar < npars; ipar++) {
    par[ipar] = f->GetParameter(ipar);
    pare[ipar] = f->GetParError(ipar);
  }

  Double_t x;
  for (Int_t ibin = 0; ibin < p->GetNbinsX(); ibin++) {
    for (Int_t itry = 0; itry < ntry; itry++) {
      for (Int_t ipar = 0; ipar < npars; ipar++)
	f->SetParameter(ipar, gRandom->Gaus(par[ipar], pare[ipar]));
      x = gRandom->Uniform(p->GetXaxis()->GetBinLowEdge(ibin + 1), p->GetXaxis()->GetBinUpEdge(ibin + 1));
      p->Fill(x, f->Eval(x));
    }
  }
}
Commit	Line	Data
aab90846	1	TH1 *
	2	HistoUtils_weightedmean(TH1 h1, TH1 h2)
	3	{
	4
	5	TH1 ho = (TH1 )h1->Clone("ho");
	6	ho->Reset();
	7	Double_t val1, val2, w1, w2, mean, meane;
	8	for (Int_t i = 0; i < ho->GetNbinsX(); i++) {
	9	val1 = h1->GetBinContent(i + 1);
	10	w1 = 1. / (h1->GetBinError(i + 1) * h1->GetBinError(i + 1));
	11	val2 = h2->GetBinContent(i + 1);
	12	w2 = 1. / (h2->GetBinError(i + 1) * h2->GetBinError(i + 1));
	13
	14	if (val1 == 0 && val2 == 0) continue;
	15
	16	mean = (w1 * val1 + w2 * val2) / (w1 + w2);
	17	meane = TMath::Sqrt(1. / (w1 + w2));
	18
	19	ho->SetBinContent(i + 1, mean);
	20	ho->SetBinError(i + 1, meane);
	21	}
	22
	23	return ho;
	24	}
	25
	26	//__________________________________________________________________
	27
	28	TH1 *
	29	HistoUtils_smartdifference(TH1 hnum, TH1 hden)
	30	{
	31
	32	TH1 hr = (TH1 )hnum->Clone("hr");
	33	hr->Reset();
	34	Double_t ref;
	35	for (Int_t i = 0; i < hr->GetNbinsX(); i++) {
	36	if (hnum->GetBinError(i + 1) <= 0.) continue;
	37	ref = hden->Interpolate(hr->GetBinCenter(i + 1));
	38	if (ref <= 0.) continue;
	39	hr->SetBinContent(i + 1, (hnum->GetBinContent(i + 1) - ref) / ref);
	40	hr->SetBinError(i + 1, hnum->GetBinError(i + 1) / ref);
	41	}
	42	return hr;
	43	}
	44
	45	//__________________________________________________________________
	46
	47	TH1 *
	48	HistoUtils_smartratio(TH1 hnum, TH1 hden)
	49	{
	50
	51	TH1 hr = (TH1 )hnum->Clone("hr");
	52	hr->Reset();
	53	Double_t ref;
	54	for (Int_t i = 0; i < hr->GetNbinsX(); i++) {
	55	if (hnum->GetBinError(i + 1) <= 0.) continue;
	56	ref = hden->Interpolate(hr->GetBinCenter(i + 1));
	57	if (ref <= 0.) continue;
	58	hr->SetBinContent(i + 1, hnum->GetBinContent(i + 1) / ref);
	59	hr->SetBinError(i + 1, hnum->GetBinError(i + 1) / ref);
	60	}
	61	return hr;
	62	}
	63
	64	//__________________________________________________________________
65
66	TH1 *
67	HistoUtils_smartratio(TH1 hnum, TGraph hden)
68	{
69
70	TH1 hr = (TH1 )hnum->Clone("hr");
71	hr->Reset();
72	Double_t ref;
73	for (Int_t i = 0; i < hr->GetNbinsX(); i++) {
74	if (hnum->GetBinError(i + 1) <= 0.) continue;
75	ref = hden->Eval(hr->GetBinCenter(i + 1));
76	if (ref <= 0.) continue;
77	hr->SetBinContent(i + 1, hnum->GetBinContent(i + 1) / ref);
78	hr->SetBinError(i + 1, hnum->GetBinError(i + 1) / ref);
79	}
80	return hr;
81	}
82
48a73d39	83	//__________________________________________________________________
	84
	85	HistoUtils_drawthemall(const Char_t *filename, Int_t sleepms = 100)
	86	{
	87	TFile *f1 = TFile::Open(filename);
	88	TList *l1 = f1->GetListOfKeys();
	89	TObject *o;
	90	Char_t *name;
	91	for (Int_t i = 0; i < l1->GetEntries(); i++) {
	92	name = l1->At(i)->GetName();
	93	o = f1->Get(name);
	94	o->Draw();
	95	gPad->Update();
	96	gSystem->Sleep(sleepms);
	97	}
	98	f1->Close();
	99	}
	100
	101	//__________________________________________________________________
	102
	103	HistoUtils_autoratio(const Char_t f1name, const Char_t f2name, const Char_t outname, const Char_t title = NULL)
	104	{
	105
	106	TFile *f1 = TFile::Open(f1name);
	107	TFile *f2 = TFile::Open(f2name);
	108	TFile *fo = TFile::Open(outname, "RECREATE");
	109	TList *l1 = f1->GetListOfKeys();
	110
	111
	112	TH1D h1, h2, *hr;
	113	Char_t *name;
	114	for (Int_t i = 0; i < l1->GetEntries(); i++) {
	115	name = l1->At(i)->GetName();
	116	h1 = (TH1D *)f1->Get(name);
	117	h2 = (TH1D *)f2->Get(name);
	118	if (!h1 \|\| !h2) continue;
	119	hr = new TH1D(*h1);
	120	hr->Divide(h2);
	121	if (title)
	122	hr->SetTitle(title);
	123	fo->cd();
	124	hr->Write();
	125	}
	126
	127	delete hr;
	128	f1->Close();
	129	f2->Close();
	130	fo->Close();
	131
	132	}
	133
	134	//__________________________________________________________________
	135
	136	HistoUtils_autosystematics(const Char_t f1name, const Char_t f2name, const Char_t outname, const Char_t title = NULL)
	137	{
	138
	139	TFile *f1 = TFile::Open(f1name);
	140	TFile *f2 = TFile::Open(f2name);
	141	if (!f1 \|\| !f1->IsOpen() \|\| !f2 \|\| !f2->IsOpen()) return;
	142	TFile *fo = TFile::Open(outname, "RECREATE");
	143	TList *l1 = f1->GetListOfKeys();
	144
	145
	146	TH1D h1, h2, *hd;
147	Char_t *name;
148	for (Int_t i = 0; i < l1->GetEntries(); i++) {
149	name = l1->At(i)->GetName();
150	h1 = (TH1D *)f1->Get(name);
151	h2 = (TH1D *)f2->Get(name);
152	if (!h1 \|\| !h2) continue;
153
154	hd = new TH1D(*h1);
155	Double_t val1, val2, vald, vale1, vale2, valde;
156	for (Int_t ii = 0; ii < hd->GetNbinsX(); ii++) {
157	val1 = h1->GetBinContent(ii + 1);
158	vale1 = h1->GetBinError(ii + 1);
159	val2 = h2->GetBinContent(ii + 1);
160	vale2 = h2->GetBinError(ii + 1);
161	if (val2 == 0.) continue;
162	vald = (val1 - val2) / val2;
163	valde = TMath::Sqrt(TMath::Abs((vale1 * vale1 - vale2 * vale2))) / val2;
164	hd->SetBinContent(ii + 1, vald);
165	hd->SetBinError(ii + 1, valde);
166	}
167
168	if (title)
169	hd->SetTitle(title);
170	fo->cd();
171	hd->Write();
172	delete hd;
173	}
174
175	f1->Close();
176	f2->Close();
177	fo->Close();
178
179	}
180
181	//__________________________________________________________________
182
183	TH1 *
184	HistoUtils_ratio(const Char_t f1name, const Char_t f2name, const Char_t h1name, const Char_t h2name)
185	{
186
187	if (!f2name) f2name = f1name;
188	TFile *f1 = TFile::Open(f1name);
189	TFile *f2 = TFile::Open(f2name);
190
191	if (!h2name) h2name = h1name;
192	TH1 h1 = (TH1 )f1->Get(h1name);
193	TH1 h2 = (TH1 )f2->Get(h2name);
194
195	TH1 *hr = h1->Clone("hr");
196	hr->Sumw2();
197	hr->Divide(h2);
198
199	return hr;
200
201	}
202
203	//__________________________________________________________________
204
205	TH1 *
206	HistoUtils_systematics(const Char_t f1name, const Char_t f2name, const Char_t h1name, const Char_t h2name)
207	{
208
209	if (!f2name) f2name = f1name;
210	TFile *f1 = TFile::Open(f1name);
211	TFile *f2 = TFile::Open(f2name);
212
213	if (!h2name) h2name = h1name;
214	TH1 h1 = (TH1 )f1->Get(h1name);
215	TH1 h2 = (TH1 )f2->Get(h2name);
216
217	TH1 *hd = h1->Clone("hd");
218	Double_t val1, val2, vald, vale1, vale2, valde;
219	for (Int_t i = 0; i < h1->GetNbinsX(); i++) {
220	val1 = h1->GetBinContent(i + 1);
221	vale1 = h1->GetBinError(i + 1);
222	val2 = h2->GetBinContent(i + 1);
223	vale2 = h2->GetBinError(i + 1);
224	if (val2 == 0.) continue;
225	vald = (val1 - val2) / val2;
226	valde = TMath::Sqrt(TMath::Abs((vale1 * vale1 - vale2 * vale2))) / val2;
227	hd->SetBinContent(i + 1, vald);
228	hd->SetBinError(i + 1, valde);
229	}
230
231	return hd;
232
233	}
234
235	//__________________________________________________________________
236
237	HistoUtils_BinLogX(TH1 *h)
238	{
239	TAxis *axis = h->GetXaxis();
240	Int_t bins = axis->GetNbins();
241	Axis_t from = axis->GetXmin();
242	Axis_t to = axis->GetXmax();
243	Axis_t width = (to - from) / bins;
244	Axis_t *new_bins = new Axis_t[bins + 1];
245	for (int i = 0; i <= bins; i++) new_bins[i] = TMath::Power(10, from + i * width);
246	axis->Set(bins, new_bins);
247	delete new_bins;
248	}
249
250	//__________________________________________________________________
251
252	HistoUtils_BinNormX(TH1 *h)
253	{
254	TAxis *axis = h->GetXaxis();
255	Int_t bins = axis->GetNbins();
256	Double_t c, ec;
257	Double_t w;
258	for (Int_t i = 0; i < bins; i++) {
259	c = h->GetBinContent(i + 1);
260	ec = h->GetBinError(i + 1);
261	w = axis->GetBinWidth(i + 1);
262	c /= w;
263	ec /= w;
264	h->SetBinContent(i + 1, c);
265	h->SetBinError(i + 1, ec);
266	}
267	}
268
269	//__________________________________________________________________
270
271	TObjArray *
272	HistoUtils_FitPeak(TF1 fitFunc, TH2 h, Float_t startSigma, Float_t nSigmaMin, Float_t nSigmaMax, Int_t minIntegral = 100., const Char_t *basename = "hParam", Bool_t monitor = kFALSE)
273	{
274
275	/* gaus function if not specified */
276	if (!fitFunc)
277	fitFunc = (TF1*)gROOT->GetFunction("gaus");
278
279	/* prepare output histos */
280	TObjArray *outArray = new TObjArray();
281	Int_t npars = fitFunc->GetNpar();
282	TH1D *hpx = h->ProjectionX("hpx");
283	TH1D *hParam;
284	for (Int_t ipar = 0; ipar < npars; ipar++) {
285	hParam = new TH1D(*hpx);
286	hParam->SetName(Form("%s_%d", basename, ipar));
287	hParam->Reset();
288	outArray->Add(hParam);
289	}
290
291	/* loop over x-bins */
292	for (Int_t ibin = 0; ibin < hpx->GetNbinsX(); ibin++) {
293
294	/* check integral */
295	if (hpx->GetBinContent(ibin + 1) < minIntegral) continue;
296	/* projection y */
297	TH1D *hpy = h->ProjectionY("hpy", ibin + 1, ibin + 1);
298	/* fit peak */
299	if (HistoUtils_FitPeak(fitFunc, hpy, startSigma, nSigmaMin, nSigmaMax) != 0) {
300	delete hpy;
301	continue;
302	}
303	/* setup output histos */
304	for (Int_t ipar = 0; ipar < npars; ipar++) {
305	hParam = (TH1D *)outArray->At(ipar);
306	hParam->SetBinContent(ibin + 1, fitFunc->GetParameter(ipar));
307	hParam->SetBinError(ibin + 1, fitFunc->GetParError(ipar));
308	}
309	/* monitor */
310	if (monitor) {
311	hpy->SetMarkerStyle(20);
312	hpy->SetMarkerColor(4);
313	hpy->Draw("E1");
314	fitFunc->Draw("same");
315	gPad->Update();
316	getchar();
317	}
318	/* delete */
319	delete hpy;
320
321	}
322	/* delete */
323	delete hpx;
324	/* return output array */
325	return outArray;
326	}
327
328	//__________________________________________________________________
329
330	Int_t
331	HistoUtils_FitPeak(TF1 fitFunc, TH1 h, Float_t startSigma, Float_t nSigmaMin, Float_t nSigmaMax)
332	{
333
334	Double_t fitCent = h->GetBinCenter(h->GetMaximumBin());
335	Double_t fitMin = fitCent - nSigmaMin * startSigma;
336	Double_t fitMax = fitCent + nSigmaMax * startSigma;
337	if (fitMin < h->GetXaxis()->GetXmin()) fitMin = h->GetXaxis()->GetXmin();
338	if (fitMax > h->GetXaxis()->GetXmax()) fitMax = h->GetXaxis()->GetXmax();
339	fitFunc->SetParameter(1, fitCent);
340	fitFunc->SetParameter(2, startSigma);
341	Int_t fitres = h->Fit(fitFunc, "WWq0", "", fitMin, fitMax);
342	if (fitres != 0) return fitres;
343	/* refit with better range */
344	for (Int_t i = 0; i < 3; i++) {
345	fitCent = fitFunc->GetParameter(1);
346	fitMin = fitCent - nSigmaMin * fitFunc->GetParameter(2);
347	fitMax = fitCent + nSigmaMax * fitFunc->GetParameter(2);
348	if (fitMin < h->GetXaxis()->GetXmin()) fitMin = h->GetXaxis()->GetXmin();
349	if (fitMax > h->GetXaxis()->GetXmax()) fitMax = h->GetXaxis()->GetXmax();
350	fitres = h->Fit(fitFunc, "q0", "", fitMin, fitMax);
351	if (fitres != 0) return fitres;
352	}
353	return fitres;
354
355	}
356
357	//__________________________________________________________________
358
359	void
360	HistoUtils_Function2Profile(TF1 fin, TProfile p, Int_t ntry = 10000)
361	{
362
363	TF1 f = new TF1(fin);
364	Int_t npars = f->GetNpar();
365	Double_t par[1000];
366	Double_t pare[1000];
367	for (Int_t ipar = 0; ipar < npars; ipar++) {
368	par[ipar] = f->GetParameter(ipar);
369	pare[ipar] = f->GetParError(ipar);
370	}
371
372	Double_t x;
373	for (Int_t ibin = 0; ibin < p->GetNbinsX(); ibin++) {
374	for (Int_t itry = 0; itry < ntry; itry++) {
375	for (Int_t ipar = 0; ipar < npars; ipar++)
376	f->SetParameter(ipar, gRandom->Gaus(par[ipar], pare[ipar]));
377	x = gRandom->Uniform(p->GetXaxis()->GetBinLowEdge(ibin + 1), p->GetXaxis()->GetBinUpEdge(ibin + 1));
378	p->Fill(x, f->Eval(x));
379	}
380	}
381	}