00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include "BAT/BCH2D.h"
00011
00012 #include "BAT/BCMath.h"
00013 #include "BAT/BCLog.h"
00014
00015 #include <TH1D.h>
00016 #include <TH2D.h>
00017 #include <TGraph.h>
00018 #include <TCanvas.h>
00019 #include <TLine.h>
00020 #include <TMarker.h>
00021 #include <TLegend.h>
00022 #include <TString.h>
00023
00024 #include <math.h>
00025
00026
00027
00028 BCH2D::BCH2D()
00029 {
00030 fHistogram = 0;
00031 fIntegratedHistogram = 0;
00032
00033 fModeFlag = 0;
00034 }
00035
00036
00037
00038 BCH2D::BCH2D(TH2D * h)
00039 {
00040 fHistogram = h;
00041 fIntegratedHistogram = 0;
00042
00043 fModeFlag = 0;
00044 }
00045
00046
00047
00048 BCH2D::~BCH2D()
00049 {
00050 if (fHistogram) delete fHistogram;
00051 if (fIntegratedHistogram) delete fIntegratedHistogram;
00052 }
00053
00054
00055
00056 void BCH2D::GetMode(double& mode)
00057 {
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067 }
00068
00069
00070
00071 void BCH2D::Print(const char * filename, int options, int ww, int wh)
00072 {
00073
00074 TCanvas * c;
00075 if(ww >0 && wh > 0)
00076 c = new TCanvas("c","c",ww,wh);
00077 else
00078 c = new TCanvas("c","c",700,700);
00079 c -> cd();
00080
00081
00082 this->Draw(options);
00083
00084 gPad->RedrawAxis();
00085
00086
00087 c -> Print(filename);
00088 }
00089
00090
00091
00092 void BCH2D::Draw(int options, bool drawmode)
00093 {
00094
00095 fHistogram -> SetLineColor(kBlack);
00096
00097
00098 fHistogram->Scale(1./fHistogram->Integral("width"));
00099
00100 double modex,modey;
00101
00102
00103 if(fModeFlag)
00104 {
00105 modex = fMode[0];
00106 modey = fMode[1];
00107 }
00108 else
00109 {
00110 int maximumbin = fHistogram -> GetMaximumBin();
00111
00112 int binx = maximumbin % (fHistogram -> GetNbinsX() + 2);
00113 int biny = maximumbin / (fHistogram -> GetNbinsX() + 2);
00114
00115 modex = fHistogram -> GetXaxis() -> GetBinCenter(binx);
00116 modey = fHistogram -> GetYaxis() -> GetBinCenter(biny);
00117 }
00118
00119
00120 fHistogram->Scale(1./fHistogram->Integral("width"));
00121
00122
00123 if (options == 0)
00124 fHistogram -> Draw("CONT0");
00125 else if (options == 1)
00126 {
00127 fHistogram -> Draw("CONT3");
00128
00129
00130 this -> CalculateIntegratedHistogram();
00131
00132 double levels[4];
00133 levels[0] = 0.;
00134 levels[1] = this -> GetLevel(1.0 - 0.6827);
00135 levels[2] = this -> GetLevel(1.0 - 0.9545);
00136 levels[3] = this -> GetLevel(1.0 - 0.9973);
00137
00138 fHistogram -> SetContour(4, levels);
00139
00140
00141 TMarker* marker = new TMarker(modex, modey, 24);
00142 marker -> Draw();
00143
00144 TLegend* legend = new TLegend(0.65, 0.80, 0.95, 0.95);
00145 legend -> SetBorderSize(0);
00146 legend -> SetFillColor(kWhite);
00147 legend -> AddEntry(fHistogram, "68% prob. region", "L");
00148 legend -> AddEntry(marker, "Best fit", "P");
00149 legend -> Draw();
00150 }
00151 else if (options == 2)
00152 {
00153 fHistogram -> Draw("CONT3");
00154
00155
00156 this -> CalculateIntegratedHistogram();
00157
00158 double levels[2];
00159 double level32 = this -> GetLevel(0.32);
00160 levels[0] = 0.;
00161 levels[1] = level32;
00162
00163 fHistogram -> SetContour(2, levels);
00164
00165
00166 TMarker* marker = new TMarker(modex, modey, 24);
00167 marker -> Draw();
00168
00169 TLegend* legend = new TLegend(0.65, 0.80, 0.95, 0.95);
00170 legend -> SetBorderSize(0);
00171 legend -> SetFillColor(kWhite);
00172 legend -> AddEntry(fHistogram, "68% prob. region", "L");
00173 legend -> AddEntry(marker, "Best fit", "P");
00174 legend -> Draw();
00175
00176 }
00177 else if (options == 3)
00178 {
00179 fHistogram -> Draw("CONT3");
00180
00181
00182 this -> CalculateIntegratedHistogram();
00183
00184 double levels[2];
00185 double level10 = this -> GetLevel(0.10);
00186 levels[0] = 0.;
00187 levels[1] = level10;
00188
00189 fHistogram -> SetContour(2, levels);
00190
00191 TLegend* legend = new TLegend(0.65, 0.80, 0.95, 0.95);
00192 legend -> SetBorderSize(0);
00193 legend -> SetFillColor(kWhite);
00194 legend -> AddEntry(fHistogram, "90% prob. region", "L");
00195 legend -> Draw();
00196 }
00197 else if (options == 4)
00198 {
00199 fHistogram -> Draw("CONT3");
00200
00201
00202 this -> CalculateIntegratedHistogram();
00203
00204 double levels[2];
00205 double level5 = this -> GetLevel(0.05);
00206 levels[0] = 0.;
00207 levels[1] = level5;
00208
00209 fHistogram -> SetContour(2, levels);
00210
00211 TLegend* legend = new TLegend(0.65, 0.80, 0.95, 0.95);
00212 legend -> SetBorderSize(0);
00213 legend -> SetFillColor(kWhite);
00214 legend -> AddEntry(fHistogram, "95% prob. region", "L");
00215 legend -> Draw();
00216 }
00217 else if (options == 5)
00218 fHistogram -> Draw("COL");
00219 else if (options == 52 || options == 521)
00220 {
00221
00222 int nx = fHistogram -> GetNbinsX();
00223 int ny = fHistogram -> GetNbinsY();
00224 double xmin = fHistogram->GetXaxis()->GetXmin();
00225 double xmax = fHistogram->GetXaxis()->GetXmax();
00226 TH2D * h = new TH2D(
00227 TString::Format("htemp52_%d",BCLog::GetHIndex()),fHistogram->GetTitle(),
00228 nx,xmin,xmax,
00229 ny,fHistogram->GetYaxis()->GetXmin(),fHistogram->GetYaxis()->GetXmax());
00230 h->SetXTitle(fHistogram->GetXaxis()->GetTitle());
00231 h->SetYTitle(fHistogram->GetYaxis()->GetTitle());
00232
00233
00234
00235 for(int i=1;i<=nx;i++)
00236 for(int j=1;j<=ny;j++)
00237 {
00238 double val = fHistogram -> GetBinContent(i,j);
00239
00240 if(options == 521)
00241 {
00242
00243
00244
00245 val = log10(val);
00246 }
00247 h -> SetBinContent(i,j,val);
00248 }
00249
00250
00251 h->SetStats(0);
00252 h->Draw("col");
00253
00254
00255 fHistogram -> Draw("cont3 same");
00256 fHistogram -> SetLineWidth(2);
00257
00258
00259 this -> CalculateIntegratedHistogram();
00260
00261 double levels[] = { this -> GetLevel(0.32) };
00262 fHistogram -> SetContour(1, levels);
00263
00264
00265 if(drawmode)
00266 {
00267 TMarker * marker0 = new TMarker(modex, modey, 8);
00268 marker0 -> SetMarkerColor(0);
00269 marker0 -> SetMarkerSize(.7);
00270 marker0 -> Draw();
00271 TMarker * marker1 = new TMarker(modex, modey, 4);
00272 marker1 -> SetMarkerColor(1);
00273 marker1 -> SetMarkerSize(.7);
00274 marker1 -> Draw();
00275
00276
00277
00278 }
00279 }
00280 else if (options == 53 || options == 531)
00281 {
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315 gPad -> SetLogz();
00316 fHistogram -> Draw("colz");
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329 TMarker * marker0 = new TMarker(modex, modey, 8);
00330 marker0 -> SetMarkerColor(0);
00331 marker0 -> SetMarkerSize(.7);
00332 marker0 -> Draw();
00333 TMarker * marker1 = new TMarker(modex, modey, 4);
00334 marker1 -> SetMarkerColor(1);
00335 marker1 -> SetMarkerSize(.7);
00336 marker1 -> Draw();
00337
00338
00339
00340 }
00341 }
00342
00343
00344
00345 void BCH2D::CalculateIntegratedHistogram()
00346 {
00347 int nz = 100;
00348
00349 double zmax = fHistogram -> GetMaximum();
00350 double dz = zmax / double(nz);
00351
00352 double nx = fHistogram -> GetNbinsX();
00353 double ny = fHistogram -> GetNbinsY();
00354
00355
00356 if (fIntegratedHistogram)
00357 delete fIntegratedHistogram;
00358
00359 fIntegratedHistogram = new TH1D(
00360 TString::Format("%s_int_prob_%d",fHistogram->GetName(),BCLog::GetHIndex()), "", nz, 0.0, zmax);
00361 fIntegratedHistogram -> SetXTitle("z");
00362 fIntegratedHistogram -> SetYTitle("Integrated probability");
00363 fIntegratedHistogram -> SetStats(kFALSE);
00364
00365
00366 for (int ix = 1; ix <= nx; ix++)
00367 for (int iy = 1; iy <= ny; iy++)
00368 {
00369 int binmin = BCMath::Nint(fHistogram -> GetBinContent(ix, iy) / dz);
00370 for (int i = binmin; i <= nz; i++)
00371 fIntegratedHistogram -> SetBinContent(i,
00372 fIntegratedHistogram -> GetBinContent(i) +
00373 fHistogram -> GetBinContent(ix, iy));
00374 }
00375 }
00376
00377
00378
00379 double BCH2D::GetLevel(double p)
00380 {
00381 double quantiles[1];
00382 double probsum[1];
00383 probsum[0] = p;
00384
00385 fIntegratedHistogram -> GetQuantiles( 1, quantiles, probsum);
00386
00387 return quantiles[0];
00388 }
00389
00390
00391
00392 TGraph ** BCH2D::GetBandGraphs(TH2D * h, int &n)
00393 {
00394 n=0;
00395
00396 int nbands=0;
00397 TH2D * hcopy = (TH2D*)h->Clone(TString::Format("%s_copy_%d",h->GetName(),BCLog::GetHIndex()));
00398
00399 std::vector <int> nint=GetNIntervalsY(hcopy,nbands);
00400
00401 if(nbands>2)
00402 {
00403 BCLog::Out(BCLog::warning,BCLog::warning,
00404 Form("BCH2D::GetBandGraphs : Detected %d bands. Maximum allowed is 2 (sorry).",nbands));
00405 return 0;
00406 }
00407 else if(nbands==0)
00408 {
00409 BCLog::Out(BCLog::warning,BCLog::warning,"BCH2D::GetBandGraphs : No bands detected.");
00410 return 0;
00411 }
00412
00413 TGraph ** gxx = new TGraph*[nbands];
00414
00415 TH2D * h0 = (TH2D*)h->Clone();
00416
00417 if (nbands>0)
00418 gxx[0] = GetLowestBandGraph(h0,nint);
00419
00420 if (nbands==2)
00421 {
00422 gxx[1] = GetLowestBandGraph(h0);
00423 n=2;
00424 }
00425 else
00426 n=1;
00427
00428 return gxx;
00429 }
00430
00431
00432
00433 std::vector <int> BCH2D::GetNIntervalsY(TH2D * h, int &nfoundmax)
00434 {
00435 std::vector <int> nint;
00436
00437 int nx = h -> GetNbinsX();
00438 int ny = h -> GetNbinsY();
00439
00440 nfoundmax=0;
00441
00442
00443 for (int ix=1; ix<=nx; ix++)
00444 {
00445 int nfound=0;
00446
00447
00448
00449 for (int iy=1; iy<=ny; iy++)
00450 if(h->GetBinContent(ix,iy)>0.)
00451 {
00452 while(h->GetBinContent(ix,++iy)>0.)
00453 ;
00454 nfound++;
00455 }
00456
00457
00458 if(nfound>nfoundmax)
00459 nfoundmax=nfound;
00460
00461 nint.push_back(nfound);
00462 }
00463
00464 return nint;
00465 }
00466
00467
00468
00469 TGraph * BCH2D::GetLowestBandGraph(TH2D * h)
00470 {
00471 int n;
00472 return GetLowestBandGraph(h,GetNIntervalsY(h,n));
00473 }
00474
00475
00476
00477 TGraph * BCH2D::GetLowestBandGraph(TH2D * h, std::vector<int> nint)
00478 {
00479 int nx = h -> GetNbinsX() - 1;
00480 int ny = h -> GetNbinsY();
00481
00482 TGraph * g = new TGraph(2*nx);
00483
00484 for (int ix=1; ix<=nx; ix++)
00485 {
00486
00487 double x;
00488 if(ix==1)
00489 x = h->GetXaxis()->GetBinLowEdge(1);
00490 else if(ix==nx)
00491 x = h->GetXaxis()->GetBinLowEdge(nx+1);
00492 else
00493 x = h->GetXaxis()->GetBinCenter(ix);
00494
00495 for(int iy=1; iy<=ny; iy++)
00496 if(h->GetBinContent(ix,iy)>0.)
00497 {
00498
00499 g->SetPoint(ix-1, x, h->GetYaxis()->GetBinLowEdge(iy));
00500
00501
00502 if(nint[ix-1]==2)
00503 h->SetBinContent(ix,iy,0.);
00504
00505 while(h->GetBinContent(ix,++iy)>0.)
00506 if(nint[ix-1]==2)
00507 h->SetBinContent(ix,iy,0.);
00508
00509
00510 g->SetPoint(2*nx-ix, x, h->GetYaxis()->GetBinLowEdge(iy));
00511
00512 break;
00513 }
00514 }
00515
00516 return g;
00517 }
00518
00519
00520
00521 std::vector <double> BCH2D::GetLevelBoundary(double level)
00522 {
00523 return GetLevelBoundary(fHistogram, level);
00524 }
00525
00526
00527
00528 std::vector <double> BCH2D::GetLevelBoundary(TH2D * h, double level)
00529 {
00530 std::vector <double> b;
00531
00532 int nx = h -> GetNbinsX();
00533
00534 b.assign(nx - 1, 0.0);
00535
00536
00537 for (int ix = 1; ix < nx; ix++)
00538 {
00539 TH1D * h1 = h -> ProjectionY(TString::Format("temphist_%d",BCLog::GetHIndex()), ix, ix);
00540
00541 int nprobSum = 1;
00542 double q[1];
00543 double probSum[] = { level };
00544
00545 h1 -> GetQuantiles(nprobSum, q, probSum);
00546
00547 b[ix-1] = q[0];
00548 }
00549
00550 return b;
00551 }
00552
00553
00554
00555 TGraph * BCH2D::GetBandGraph(double l1, double l2)
00556 {
00557 return GetBandGraph(fHistogram , l1, l2);
00558 }
00559
00560
00561
00562 TGraph * BCH2D::GetBandGraph(TH2D * h , double l1, double l2)
00563 {
00564
00565 int nx = h -> GetNbinsX() - 1;
00566
00567 TGraph * g = new TGraph(2*nx);
00568
00569
00570
00571
00572 std::vector <double> ymin = GetLevelBoundary(h,l1);
00573 std::vector <double> ymax = GetLevelBoundary(h,l2);
00574
00575 for (int i = 0; i < nx; i++)
00576 {
00577 g -> SetPoint(i, h -> GetXaxis() -> GetBinCenter(i+1), ymin[i]);
00578 g -> SetPoint(nx+i, h -> GetXaxis() -> GetBinCenter(nx-i), ymax[nx-i-1]);
00579 }
00580
00581 return g;
00582 }
00583
00584