PWGLF/SPECTRA/PiKaPr/TOF/pPb502/macros/HistoUtils.C

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