PWGLF/SPECTRA/PiKaPr/TOF/PbPb276/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, 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
  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 }