| 1 |
#include <iostream>
|
| 2 |
#include <iomanip>
|
| 3 |
#include <TRandom.h>
|
| 4 |
#include <TF1.h>
|
| 5 |
#include <TH1.h>
|
| 6 |
#include <TCanvas.h>
|
| 7 |
#include <TMath.h>
|
| 8 |
#include <sstream>
|
| 9 |
using namespace std;
|
| 10 |
|
| 11 |
void find68(TH1F* histo, long N, int rounds)
|
| 12 |
{
|
| 13 |
float max = histo->GetBinCenter(histo->GetMaximumBin());
|
| 14 |
int startbin=histo->GetMaximumBin();
|
| 15 |
int left=startbin;
|
| 16 |
int right=startbin;
|
| 17 |
float want=(float)N * 0.68*rounds;
|
| 18 |
float have=histo->GetBinContent(startbin);
|
| 19 |
cout << "found maximum at " << max << " ( this is bin " << startbin << ")" << endl;
|
| 20 |
while(have<want)
|
| 21 |
{
|
| 22 |
if(histo->GetBinContent(left-1)>histo->GetBinContent(right+1))
|
| 23 |
{
|
| 24 |
have+=histo->GetBinContent(left-1);
|
| 25 |
left--;
|
| 26 |
}
|
| 27 |
else
|
| 28 |
{
|
| 29 |
have+=histo->GetBinContent(right+1);
|
| 30 |
right++;
|
| 31 |
}
|
| 32 |
}
|
| 33 |
float loweredge=histo->GetBinCenter(left);
|
| 34 |
float upperedge=histo->GetBinCenter(right);
|
| 35 |
cout << "we now have " << have << " and we wanted at least " << want << endl;
|
| 36 |
cout << "the range is " << loweredge << " to " << upperedge<< " corresponding to :" << endl;
|
| 37 |
|
| 38 |
cout << "\033[1;34m " << max << "+ " << (upperedge-max) << " - " << (max-loweredge) << " \033[0m" << endl;
|
| 39 |
|
| 40 |
// cout << max << "+ " << (upperedge-max) << " - " << (max-loweredge) << endl;
|
| 41 |
}
|
| 42 |
|
| 43 |
|
| 44 |
void combine_dists(TH1F *combineddist,float *x1,float *x2,float *x3, long N)
|
| 45 |
{
|
| 46 |
float temp=-99;
|
| 47 |
for (int i=0;i<N;i++)
|
| 48 |
{
|
| 49 |
temp = x1[i]+x2[i]-x3[i];
|
| 50 |
combineddist->Fill(temp);
|
| 51 |
}
|
| 52 |
}
|
| 53 |
|
| 54 |
void generate(float *x1, int lambda,long N,TH1F* histo)
|
| 55 |
{
|
| 56 |
stringstream function;
|
| 57 |
function << "(TMath::Power(" << lambda << ",x)*TMath::Exp(-" << lambda << "))/(TMath::Gamma(x+1))";
|
| 58 |
//TF1 *f1 = new TF1("f1","(TMath::Power(1,x)*TMath::Exp(-1))/(TMath::Gamma(x+1))",0,10);
|
| 59 |
TF1 *f1 = new TF1("f1",function.str().c_str(),0,10);
|
| 60 |
for(int i=0;i<N;i++)
|
| 61 |
{
|
| 62 |
x1[i]=f1->GetRandom();
|
| 63 |
histo->Fill(x1[i]);
|
| 64 |
}//end of for
|
| 65 |
}//end of generate_x1
|
| 66 |
|
| 67 |
|
| 68 |
int main()
|
| 69 |
{
|
| 70 |
long N=1000000;
|
| 71 |
int rounds=1;
|
| 72 |
cout << "generating the distributions with N= " << N << endl;
|
| 73 |
TH1F *histo = new TH1F("histo","histo",100,-10,10);
|
| 74 |
float x1[N];
|
| 75 |
float x2[N];
|
| 76 |
float x3[N];
|
| 77 |
|
| 78 |
TH1F *summary = new TH1F("summary","summary",100,0,20);
|
| 79 |
TCanvas *c1 = new TCanvas("c1","c1");
|
| 80 |
|
| 81 |
for (int i=0;i<rounds;i++)
|
| 82 |
{
|
| 83 |
TH1F *combineddist = new TH1F("combineddist","combineddist",100,0,20);
|
| 84 |
generate(x1,5,N,histo);
|
| 85 |
generate(x2,4,N,histo);
|
| 86 |
generate(x3,1,N,histo);
|
| 87 |
combine_dists(combineddist,x1,x2,x3,N);
|
| 88 |
summary->Add(combineddist);
|
| 89 |
summary->Draw();
|
| 90 |
|
| 91 |
stringstream saveasname;
|
| 92 |
saveasname << "summary_for_i_" << i << ".png";
|
| 93 |
// c1->SaveAs(saveasname.str().c_str());
|
| 94 |
|
| 95 |
delete combineddist;
|
| 96 |
cout << "Currently there are " << summary->GetEntries() << " entries" << endl;
|
| 97 |
}
|
| 98 |
|
| 99 |
// combine_dists(combineddist2,x1,x2,x3,N);
|
| 100 |
|
| 101 |
summary->Draw();
|
| 102 |
c1->SaveAs("Canvas.png");
|
| 103 |
|
| 104 |
find68(summary,N,rounds);
|
| 105 |
|
| 106 |
return 0;
|
| 107 |
}
|
| 108 |
|
