PWGLF/SPECTRA/UTILS/YieldMean.C

   1 /* definition of the fields in the histogram returned */
   2 enum EValue_t {
   3   kYield = 1,
   4   kYieldStat,
   5   kYieldSysHi,
   6   kYieldSysLo,
   7   kMean,
   8   kMeanStat,
   9   kMeanSysHi,
  10   kMeanSysLo
  11 };
  12
  13 TH1 *
  14 YieldMean(TH1 *hstat, TH1 *hsys, TF1 *f = NULL, Double_t min = 0., Double_t max = 10., Double_t loprecision = 0.01, Double_t hiprecision = 0.1, Option_t *opt = "0q",TString logfilename="log.root")
  15 {
  16   /* set many iterations when fitting the data so we don't
  17      stop minimization with MAX_CALLS */
  18   TVirtualFitter::SetMaxIterations(1000000);
  19
  20   /* create output histo */
  21   Double_t integral, mean;
  22   TH1 *hout = new TH1D("hout", "", 8, 0, 8);
  23   TH1 *hlo, *hhi;
  24
  25   /* create histo with stat+sys errors */
  26   TH1 *htot = (TH1 *)hstat->Clone(Form("%sfittedwith%s",hstat->GetName(),f->GetName()));
  27   for (Int_t ibin = 0; ibin < htot->GetNbinsX(); ibin++) {
  28     htot->SetBinError(ibin + 1, TMath::Sqrt(hsys->GetBinError(ibin + 1) * hsys->GetBinError(ibin + 1) + hstat->GetBinError(ibin + 1) * hstat->GetBinError(ibin + 1)));
  29   }
  30
  31   /*
  32    *   measure the central value
  33    */
  34
  35   do Int_t fitres = htot->Fit(f, opt);
  36   while (fitres != 0);
  37   TFile* filewithfits=TFile::Open(logfilename.Data(),"UPDATE");
  38   htot->Write();
  39   filewithfits->Close();
  40   delete filewithfits;
  41
  42   cout<<" Fit sys+stat for " <<f->GetName()<<endl;
  43   cout<<"NDF="<<f->GetNDF()<<" Chi^2="<<f->GetChisquare()<<" Chi^2/NDF="<<f->GetChisquare()/f->GetNDF()<<endl;
  44
  45   hlo = YieldMean_LowExtrapolationHisto(htot, f, min, loprecision);
  46   hhi = YieldMean_HighExtrapolationHisto(htot, f, max, hiprecision);
  47   YieldMean_IntegralMean(htot, hlo, hhi, integral, mean,kTRUE);
  48   hout->SetBinContent(kYield, integral);
  49   hout->SetBinContent(kMean, mean);
  50
  51   /*
  52    * STATISTICS
  53    */
  54
  55   TCanvas *cCanvasStat = new TCanvas("cCanvasStat");
  56   cCanvasStat->Divide(2, 1);
  57
  58   /*
  59    * measure statistical error
  60    */
  61
  62   /* fit with stat error */
  63   do Int_t fitres = hstat->Fit(f, opt);
  64   while (fitres != 0);
  65   hlo = YieldMean_LowExtrapolationHisto(hstat, f, min, loprecision);
  66   hhi = YieldMean_HighExtrapolationHisto(hstat, f, max, hiprecision);
  67
  68   /* random generation with integration (coarse) */
  69   TH1 *hIntegral_tmp = new TH1F("hIntegral_tmp", "", 1000, 0.75 * integral, 1.25 * integral);
  70   TH1 *hMean_tmp = new TH1F("hMean_tmp", "", 1000, 0.75 * mean, 1.25 * mean);
  71   for (Int_t irnd = 0; irnd < 100; irnd++) {
  72     /* get random histogram */
  73     TH1 *hrnd = YieldMean_ReturnRandom(hstat);
  74     /* fit */
  75     TH1 *hrndlo = YieldMean_ReturnCoherentRandom(hlo);
  76     TH1 *hrndhi = YieldMean_ReturnCoherentRandom(hhi);
  77     /* integrate */
  78     YieldMean_IntegralMean(hrnd, hrndlo, hrndhi, integral, mean);
  79     hIntegral_tmp->Fill(integral);
  80     hMean_tmp->Fill(mean);
  81     delete hrnd;
  82     delete hrndlo;
  83     delete hrndhi;
  84   }
  85   /* random generation with integration (fine) */
  86   TH1 *hIntegral = new TH1F("hIntegral", "", 100,
  87                             hIntegral_tmp->GetMean() - 10. * hIntegral_tmp->GetRMS(),
  88                             hIntegral_tmp->GetMean() + 10. * hIntegral_tmp->GetRMS());
  89   TH1 *hMean = new TH1F("hMean", "", 100,
  90                         hMean_tmp->GetMean() - 10. * hMean_tmp->GetRMS(),
  91                         hMean_tmp->GetMean() + 10. * hMean_tmp->GetRMS());
  92   for (Int_t irnd = 0; irnd < 1000; irnd++) {
  93     /* get random histogram */
  94     TH1 *hrnd = YieldMean_ReturnRandom(hstat);
  95     /* fit */
  96     TH1 *hrndlo = YieldMean_ReturnCoherentRandom(hlo);
  97     TH1 *hrndhi = YieldMean_ReturnCoherentRandom(hhi);
  98     /* integrate */
  99     YieldMean_IntegralMean(hrnd, hrndlo, hrndhi, integral, mean);
 100     hIntegral->Fill(integral);
 101     hMean->Fill(mean);
 102     delete hrnd;
 103     delete hrndlo;
 104     delete hrndhi;
 105   }
 106   TF1 *gaus = (TF1 *)gROOT->GetFunction("gaus");
 107
 108   cCanvasStat->cd(1);
 109   hIntegral->Fit(gaus, "q");
 110   integral = hout->GetBinContent(kYield) * gaus->GetParameter(2) / gaus->GetParameter(1);
 111   hout->SetBinContent(kYieldStat, integral);
 112
 113   cCanvasStat->cd(2);
 114   hMean->Fit(gaus, "q");
 115   mean = hout->GetBinContent(kMean) * gaus->GetParameter(2) / gaus->GetParameter(1);
 116   hout->SetBinContent(kMeanStat, mean);
 117
 118   /*
 119    * SYSTEMATICS
 120    */
 121
 122   TCanvas *cCanvasSys = new TCanvas("cCanvasYieldSys");
 123   cCanvasSys->Divide(2, 1);
 124   cCanvasSys->cd(1)->DrawFrame(min, 1.e-3, max, 1.e3);
 125   hsys->SetMarkerStyle(20);
 126   hsys->SetMarkerColor(1);
 127   hsys->SetMarkerSize(1);
 128   hsys->Draw("same");
 129   cCanvasSys->cd(2)->DrawFrame(min, 1.e-3, max, 1.e3);
 130   hsys->Draw("same");
 131
 132   /*
 133    * systematic error high
 134    */
 135
 136   TH1 *hhigh = YieldMean_ReturnExtremeHighHisto(hsys);
 137   do Int_t fitres = hhigh->Fit(f, opt);
 138   while (fitres != 0);
 139   hlo = YieldMean_LowExtrapolationHisto(hhigh, f, min, loprecision);
 140   hhi = YieldMean_HighExtrapolationHisto(hhigh, f, max, hiprecision);
 141   YieldMean_IntegralMean(hhigh, hlo, hhi, integral, mean);
 142   integral = TMath::Abs(integral - hout->GetBinContent(kYield));
 143   hout->SetBinContent(kYieldSysHi, integral);
 144
 145   cCanvasSys->cd(1);
 146   f->SetLineColor(2);
 147   f->DrawCopy("same");
 148
 149   /*
 150    * systematic error hard
 151    */
 152
 153   TH1 *hhard = YieldMean_ReturnExtremeHardHisto(hsys);
 154   do Int_t fitres = hhard->Fit(f, opt);
 155   while (fitres != 0);
 156   hlo = YieldMean_LowExtrapolationHisto(hhard, f, min, loprecision);
 157   hhi = YieldMean_HighExtrapolationHisto(hhard, f, max, hiprecision);
 158   YieldMean_IntegralMean(hhard, hlo, hhi, integral, mean);
 159   mean = TMath::Abs(mean - hout->GetBinContent(kMean));
 160   hout->SetBinContent(kMeanSysHi, mean);
 161
 162   cCanvasSys->cd(2);
 163   f->SetLineColor(2);
 164   f->DrawCopy("same");
 165
 166   /*
 167    * systematic error low
 168    */
 169
 170   TH1 *hlow = YieldMean_ReturnExtremeLowHisto(hsys);
 171   do Int_t fitres = hlow->Fit(f, opt);
 172   while (fitres != 0);
 173   hlo = YieldMean_LowExtrapolationHisto(hlow, f, min, loprecision);
 174   hhi = YieldMean_HighExtrapolationHisto(hlow, f, max, hiprecision);
 175   YieldMean_IntegralMean(hlow, hlo, hhi, integral, mean);
 176   integral = TMath::Abs(integral - hout->GetBinContent(kYield));
 177   hout->SetBinContent(kYieldSysLo, integral);
 178
 179   cCanvasSys->cd(1);
 180   f->SetLineColor(4);
 181   f->DrawCopy("same");
 182
 183   /*
 184    * systematic error soft
 185    */
 186
 187   TH1 *hsoft = YieldMean_ReturnExtremeSoftHisto(hsys);
 188   do Int_t fitres = hsoft->Fit(f, opt);
 189   while (fitres != 0);
 190   hlo = YieldMean_LowExtrapolationHisto(hsoft, f, min, loprecision);
 191   hhi = YieldMean_HighExtrapolationHisto(hsoft, f, max, hiprecision);
 192   YieldMean_IntegralMean(hsoft, hlo, hhi, integral, mean);
 193   mean = TMath::Abs(mean - hout->GetBinContent(kMean));
 194   hout->SetBinContent(kMeanSysLo, mean);
 195
 196   cCanvasSys->cd(2);
 197   f->SetLineColor(4);
 198   f->DrawCopy("same");
 199
 200   return hout;
 201 }
 202
 203 TH1 *
 204 YieldMean_LowExtrapolationHisto(TH1 *h, TF1 *f, Double_t min, Double_t binwidth = 0.01)
 205 {
 206   /* find lowest edge in histo */
 207   Int_t binlo;
 208   Double_t lo;
 209   for (Int_t ibin = 1; ibin < h->GetNbinsX() + 1; ibin++) {
 210     if (h->GetBinContent(ibin) != 0.) {
 211       binlo = ibin;
 212       lo = h->GetBinLowEdge(ibin);
 213       break;
 214     }
 215   }
 216
 217   Int_t nbins = (lo - min) / binwidth;
 218   TH1 *hlo = new TH1F("hlo", "", nbins, min, lo);
 219
 220   /* integrate function in histogram bins */
 221   Double_t cont, err, width;
 222   for (Int_t ibin = 0; ibin < hlo->GetNbinsX(); ibin++) {
 223     width = hlo->GetBinWidth(ibin + 1);
 224     cont = f->Integral(hlo->GetBinLowEdge(ibin + 1), hlo->GetBinLowEdge(ibin + 2), (Double_t *)0, 1.e-6);
 225     err = f->IntegralError(hlo->GetBinLowEdge(ibin + 1), hlo->GetBinLowEdge(ibin + 2), (Double_t *)0, (Double_t *)0, 1.e-6);
 226     hlo->SetBinContent(ibin + 1, cont / width);
 227     hlo->SetBinError(ibin + 1, err / width);
 228   }
 229
 230   return hlo;
 231 }
 232
 233 TH1 *
 234 YieldMean_HighExtrapolationHisto(TH1 *h, TF1 *f, Double_t max, Double_t binwidth = 0.1)
 235 {
 236   /* find highest edge in histo */
 237   Int_t binhi;
 238   Double_t hi;
 239   for (Int_t ibin = h->GetNbinsX(); ibin > 0; ibin--) {
 240     if (h->GetBinContent(ibin) != 0.) {
 241       binhi = ibin + 1;
 242       hi = h->GetBinLowEdge(ibin + 1);
 243       break;
 244     }
 245   }
 246
 247   Int_t nbins = (max - hi) / binwidth;
 248   TH1 *hhi = new TH1F("hhi", "", nbins, hi, max);
 249
 250   /* integrate function in histogram bins */
 251   Double_t cont, err, width;
 252   for (Int_t ibin = 0; ibin < hhi->GetNbinsX(); ibin++) {
 253     width = hhi->GetBinWidth(ibin + 1);
 254     cont = f->Integral(hhi->GetBinLowEdge(ibin + 1), hhi->GetBinLowEdge(ibin + 2), (Double_t *)0, 1.e-6);
 255     err = f->IntegralError(hhi->GetBinLowEdge(ibin + 1), hhi->GetBinLowEdge(ibin + 2), (Double_t *)0, (Double_t *)0, 1.e-6);
 256     hhi->SetBinContent(ibin + 1, cont / width);
 257     hhi->SetBinError(ibin + 1, err / width);
 258   }
 259
 260   return hhi;
 261 }
 262
 263 TH1 *
 264 YieldMean_ReturnRandom(TH1 *hin)
 265 {
 266   TH1 *hout = (TH1 *)hin->Clone("hout");
 267   hout->Reset();
 268   Double_t cont, err;
 269   for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
 270     if (hin->GetBinError(ibin + 1) <= 0.) continue;
 271     cont = hin->GetBinContent(ibin + 1);
 272     err = hin->GetBinError(ibin + 1);
 273     hout->SetBinContent(ibin + 1, gRandom->Gaus(cont, err));
 274     hout->SetBinError(ibin + 1, err);
 275   }
 276   return hout;
 277 }
 278
 279 TH1 *
 280 YieldMean_ReturnCoherentRandom(TH1 *hin)
 281 {
 282   TH1 *hout = (TH1 *)hin->Clone("hout");
 283   hout->Reset();
 284   Double_t cont, err, cohe;
 285   cohe = gRandom->Gaus(0., 1.);
 286   for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
 287     if (hin->GetBinError(ibin + 1) <= 0.) continue;
 288     cont = hin->GetBinContent(ibin + 1);
 289     err = hin->GetBinError(ibin + 1);
 290     hout->SetBinContent(ibin + 1, cont + cohe * err);
 291     hout->SetBinError(ibin + 1, err);
 292   }
 293   return hout;
 294 }
 295
 296 TH1 *
 297 YieldMean_ReturnExtremeHighHisto(TH1 *hin)
 298 {
 299   TH1 *hout = (TH1 *)hin->Clone(Form("%s_extremehigh", hin->GetName()));
 300   for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
 301     if (hin->GetBinError(ibin + 1) <= 0.) continue;
 302     Double_t val = hin->GetBinContent(ibin + 1);
 303     Double_t err = hin->GetBinError(ibin + 1);
 304     hout->SetBinContent(ibin + 1, val + err);
 305   }
 306   return hout;
 307 }
 308
 309 TH1 *
 310 YieldMean_ReturnExtremeLowHisto(TH1 *hin)
 311 {
 312   TH1 *hout = (TH1 *)hin->Clone(Form("%s_extremelow", hin->GetName()));
 313   for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
 314     if (hin->GetBinError(ibin + 1) <= 0.) continue;
 315     Double_t val = hin->GetBinContent(ibin + 1);
 316     Double_t err = hin->GetBinError(ibin + 1);
 317     hout->SetBinContent(ibin + 1, val - err);
 318   }
 319   return hout;
 320 }
 321
 322 TH1 *
 323 YieldMean_ReturnExtremeSoftHisto(TH1 *hin)
 324 {
 325   return YieldMean_ReturnExtremeHisto(hin, -1.);
 326 }
 327
 328 TH1 *
 329 YieldMean_ReturnExtremeHardHisto(TH1 *hin)
 330 {
 331   return YieldMean_ReturnExtremeHisto(hin, 1.);
 332 }
 333
 334 TH1 *
 335 YieldMean_ReturnExtremeHisto(TH1 *hin, Float_t sign = 1.)
 336 {
 337   Double_t ptlow, pthigh;
 338   for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
 339     if (hin->GetBinError(ibin + 1) <= 0.) continue;
 340     ptlow = hin->GetBinLowEdge(ibin + 1);
 341     break;
 342   }
 343   for (Int_t ibin = hin->GetNbinsX(); ibin >= 0; ibin--) {
 344     if (hin->GetBinError(ibin + 1) <= 0.) continue;
 345     pthigh = hin->GetBinLowEdge(ibin + 2);
 346     break;
 347   }
 348
 349   Double_t mean = hin->GetMean();
 350   Double_t maxdiff = 0.;
 351   TH1 *hmax = NULL;
 352   for (Int_t inode = 0; inode < hin->GetNbinsX(); inode++) {
 353
 354     Double_t ptnode = hin->GetBinCenter(inode + 1);
 355     TH1 *hout = (TH1 *)hin->Clone(Form("%s_extremehard", hin->GetName()));
 356
 357     for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) {
 358       if (hin->GetBinError(ibin + 1) <= 0.) continue;
 359       Double_t val = hin->GetBinContent(ibin + 1);
 360       Double_t err = hin->GetBinError(ibin + 1);
 361       Double_t cen = hin->GetBinCenter(ibin + 1);
 362       if (cen < ptnode)
 363         err *= -1. + (cen - ptlow) / (ptnode - ptlow);
 364       else
 365         err *= (cen - ptnode) / (pthigh - ptnode);
 366
 367       hout->SetBinContent(ibin + 1, val + sign * err);
 368     }
 369
 370     Double_t diff = TMath::Abs(mean - hout->GetMean());
 371     if (diff > maxdiff) {
 372       //      printf("found max at %f\n", ptnode);
 373       if (hmax) delete hmax;
 374       hmax = (TH1 *)hout->Clone("hmax");
 375       maxdiff = diff;
 376     }
 377     delete hout;
 378   }
 379   return hmax;
 380 }
 381
 382 YieldMean_IntegralMean(TH1 *hdata, TH1 *hlo, TH1 *hhi, Double_t &integral, Double_t &mean,Bool_t printinfo=kFALSE)
 383 {
 384
 385   /*
 386    * compute integrals
 387    */
 388
 389   Double_t cont, err, width, cent;
 390   Double_t I = 0., IX = 0., Ierr = 0., IXerr = 0., Ilerr = 0., IXlerr = 0.;
 391   Double_t M = 0., Merr = 0., Mlerr = 0., C;
 392   Double_t dataonly=0.0;
 393
 394   /* integrate the data */
 395   for (Int_t ibin = 0; ibin < hdata->GetNbinsX(); ibin++) {
 396     cent = hdata->GetBinCenter(ibin + 1);
 397     width = hdata->GetBinWidth(ibin + 1);
 398     cont = width * hdata->GetBinContent(ibin + 1);
 399     err = width * hdata->GetBinError(ibin + 1);
 400     if (err <= 0.) continue;
 401     I += cont;
 402     IX += cont * cent;
 403   }
 404   dataonly=I;
 405   /* integrate low */
 406   for (Int_t ibin = 0; ibin < hlo->GetNbinsX(); ibin++) {
 407     cent = hlo->GetBinCenter(ibin + 1);
 408     width = hlo->GetBinWidth(ibin + 1);
 409     cont = width * hlo->GetBinContent(ibin + 1);
 410     err = width * hlo->GetBinError(ibin + 1);
 411     if (err <= 0.) continue;
 412     I += cont;
 413     IX += cont * cent;
 414   }
 415   /* integrate high */
 416   for (Int_t ibin = 0; ibin < hhi->GetNbinsX(); ibin++) {
 417     cent = hhi->GetBinCenter(ibin + 1);
 418     width = hhi->GetBinWidth(ibin + 1);
 419     cont = width * hhi->GetBinContent(ibin + 1);
 420     err = width * hhi->GetBinError(ibin + 1);
 421     if (err <= 0.) continue;
 422     I += cont;
 423     IX += cont * cent;
 424   }
 425
 426   /* set values */
 427   integral = I;
 428   mean = IX / I;
 429   if(printinfo)
 430         cout<<"data only = "<<dataonly<<" total = "<<I<<" ratio= "<<dataonly/I<<endl;
 431 }