/*************************************************************************** * -------------------------------------------------------------------------- * Solution to exercise Random Numbers and CLT * -> Book and fill a histogram with 10000 uniform random numbers * -------------------------------------------------------------------------- * 13.06.2005, Klaus Rabbertz * * modification log: * kr, 04.07.2006: Removed explicit ex10 naming * as, 06.06.2008: Some fixes for running as macro * gq, 20.12.2009: added output of correlation ***************************************************************************/ /* Note that this file can be either used as a compiled program * or as a ROOT macro. * If it is used as a compiled program, additional include statements * and the definition of the main program have to be made. This is * not needed if the code is used as a macro. */ #ifndef __CINT__ // These include-statements are needed if the program is // run as a "stand-alone application", i.e. if it is not // called from an interactive ROOT session #include #include "TFile.h" #include "TH1.h" #include "TH2.h" #include "TProfile.h" #include "TRandom.h" #include "TMath.h" #include "TCanvas.h" #include "TPad.h" #include "TF1.h" #include "TStyle.h" #include "TText.h" #include void rancor(); void erprop(); void ex11_lsg(); //______________________________________________________________________________ int main() { ex11_lsg(); return 0; } #endif //______________________________________________________________________________ // from here on, the code can also be used as a macro void ex11_lsg() { rancor(); erprop(); } Double_t NormBinomial(Double_t *x, Double_t *par) // Binomial distribution { int n,k; double p,norm; n=par[1]; k=x[0]; p=par[0]; norm=par[2]; return norm*TMath::Binomial(n,k)*TMath::Power(p,k)*TMath::Power(1-p,n-k); } Double_t NormPoisson(Double_t *x, Double_t *par) // Poisson distribution { double norm=par[1], mu=par[0]; return norm*TMath::Poisson(x[0],mu); } void rancor() { const int Nrndm=10000; // Create a new ROOT file. TFile hfile("rancor.root","RECREATE","rancor"); // Book a 1D histogram "hran" (or any other name) holding 100 random numbers. TH1F *hran = new TH1F("hran","uniform random numbers",5, 0.0,1.0); // Create vector of pointers to 1D histograms. // Each pointer points to a histogram which holds the content of a // specific bin. TH1F *hlist[5]; // Book 1D histograms for content of each bin of the above histogram hlist[0] = new TH1F("hran01","bin 1",50,0.0,50.0); hlist[1] = new TH1F("hran02","bin 2",50,0.0,50.0); hlist[2] = new TH1F("hran03","bin 3",50,0.0,50.0); hlist[3] = new TH1F("hran04","bin 4",50,0.0,50.0); hlist[4] = new TH1F("hran05","bin 5",50,0.0,50.0); /* book a 2D histogram for the scatter plot * TH2F is a ROOT class for 2D histograms hscat in an instance of that class * Arguments: TH2F("identifier", "some title", NumberBinsX, XMin, XMax, * NumberBinsY, YMin, YMax); */ TH2F *hscat = new TH2F("hscat","scatter plot: bin 2 vs bin 4", 50, 0.0, 50.0, 50,0.0, 50.0); // Book a profile histogram // TProfile is a ROOT class for profile histograms, hprof is an instance // of that class. // Arguments: // TProfile("identifier", "some title", NumberBins, XMin, XMax, YMin, YMax, Option); // Note that YMin and YMax can be left blank/ // The Option sets the treatment of errors (refer to the ROOT manual) // "Option" can be left blank, too. TProfile *hprof = new TProfile("hprof","Profile histogram bin 2 vs bin 4",50,0.0,50.0); // Now all histograms have been booked and we can start the real work. Float_t random; for (int j=0; jReset(); // Generate 100 random numbers for each experiment. for (int i=0; i<100; i++) { random = gRandom->Rndm(); // generate uniform random number hran->Fill(random); // fill random number into histogram } // Unpack histogram into array content. Float_t *content = hran->GetArray(); // create a pointer which points to the // content of each bin in histogram h100 /* Note that the array "content" is organised as follows: * content[0] = number of underflows * content[1] = number of entries in first bin * content[2] = number of entries in second bin * .... * content[n] = number of entries last bin (assuming that the histogram has n bins) * content[n+1] = number of overflows */ for (int k=0;k<5;k++) { // loop over all bins hlist[k]->Fill(content[k+1]); // fill content of each bin into histogram } // Fill the scatter plot with the content of bin 2 and bin 4. hscat->Fill(content[2],content[4]); // Fill the profile plot with the content of bin 2 and bin 4. hprof->Fill(content[2],content[4]); } // print correlation coefficient at the end std::cout << "Correlation coefficient: "<< hscat->GetCorrelationFactor() <Draw("surf3"); //hprof->Draw(); // define Function to overlay Binomial Distribution TF1 *binom = new TF1("binom",NormBinomial,0.,50.,50); binom->SetParameter(0,0.2); binom->SetParameter(1,100.); binom->SetParameter(2,Nrndm); TCanvas c("c","Canvas",0, 0, 500, 500); hlist[0]->Draw(); binom->DrawCopy("same"); // Save all objects in the ROOT file and close it. c.Write(); hfile.Write(); hfile.Close(); } void erprop() { // Teste Fehlerfortpflanzung G. Quast, Dez. 2010 // Fülle Quotienten und Differenz Gauss-verteilter Zufallszahlen // in Histogram und vergleiche mit Gauss-Verteilungen const double mut=1.0; const double muy=0.5; const double sigt=0.3; const double sigy=0.1; const int Nrndm=10000; // eventually, create a new ROOT file. // TFile hfile("erprop.root","RECREATE","erprop"); // * define some constants const int Nbins=100; const double xmin=-0.5; const double xmax=4.5; double t,y; // ---------------------------------------------------------------- // create canvas with pads // at xtop ytop sizex sizey TCanvas *c1=new TCanvas("c1","Canvas for test of error propagation" , 100, 100, 600, 600); TPad *p1=new TPad("Pad3","Pad3",0.0,0.5,0.5,1.); p1->Draw(); TPad *p2=new TPad("Pad4","Pad4",0.5,0.5,1.,1.); p2->Draw(); TPad *p3=new TPad("Pad1","Pad1",0.,0.,0.5,0.5); p3->Draw(); TPad *p4=new TPad("Pad2","Pad2",0.5,0.0,1.,0.5); p4->Draw(); // set as active pad c1->cd(); gPad->SetFillColor(0); gPad->SetGrid(); // some histogram style options gStyle->SetHistLineColor(46); gStyle->SetHistFillColor(3); gStyle->SetHistLineWidth(3); gStyle->SetAxisColor(4,"X"); gStyle->SetAxisColor(4,"Y"); gStyle->SetLabelColor(4,"X"); gStyle->SetLabelColor(4,"Y"); // book a histogram, w. Nbins bins, range [xmin,xmax[ TH1F* htovery =new TH1F("t/y","ratio of t and y",Nbins,xmin,xmax); TH1F* hyovert =new TH1F("y/t","ratio of y and t",Nbins,xmin,xmax/2); TH1F* htminusy =new TH1F("t-y","difference ",Nbins,xmin,xmax/2); // fill random numbers in a loop , Nfill times, Nrndm at a time for(int j =0; jGaus(); y=muy+sigy*gRandom->Gaus(); htovery->Fill(t/y); hyovert->Fill(y/t); htminusy->Fill(t-y); } // normalised Gauss function TF1 *Ngauss=new TF1("gauss","[2]/sqrt(2.*pi)/[1]*exp(-(x-[0])*(x-[0])/(2.*[1]*[1]))",xmin,xmax); char txt[50]; double mu0; double sig0; TText ttxt; // draw text, histograms and Gauss functions into pads prepared above p4->cd(); ttxt.SetTextSize(0.07); sprintf(txt,"Y: mu=%1.3g sig=%1.3g",muy,sigy); ttxt.DrawText(0.1,0.8, const_cast(txt)); sprintf(txt,"T: mu=%1.3g sig=%1.3g",mut,sigt); ttxt.DrawText(0.1,0.7, const_cast(txt)); ttxt.SetTextSize(0.05); p1->cd(); htovery->Draw(); mu0=mut/muy; sig0=sqrt(sigt/mut*sigt/mut+sigy/muy*sigy/muy)*mut/muy; Ngauss->SetParameter(0, mu0); Ngauss->SetParameter(1, sig0); Ngauss->SetParameter(2, Nrndm * (xmax-xmin)/Nbins); Ngauss->DrawCopy("same"); p4->cd(); sprintf(txt,"T/Y Gauss: mu=%1.3g sig=%1.3g ==> S1=%1.2g", mu0,sig0,fabs(1.-mu0)/sig0); ttxt.DrawText(0.1,0.6, const_cast(txt)); sprintf(txt," ToyMC: mean=%1.3g RMS=%1.3g",htovery->GetMean(),htovery->GetRMS()); ttxt.DrawText(0.15,0.55, const_cast(txt)); p2->cd(); hyovert->Draw(); mu0=muy/mut; sig0=sqrt(sigt/mut*sigt/mut+sigy/muy*sigy/muy)*muy/mut; Ngauss->SetParameter(0, mu0); Ngauss->SetParameter(1, sig0); Ngauss->SetParameter(2, Nrndm * (xmax/2-xmin)/Nbins); Ngauss->DrawCopy("same"); p4->cd(); sprintf(txt,"Y/T Gauss: mu=%1.3g sig=%1.3g ==> S2=%1.2g", mu0,sig0,fabs(1.-mu0)/sig0); ttxt.DrawText(0.1,0.4, const_cast(txt)); sprintf(txt," ToyMC: mean=%1.3g RMS=%1.3g",hyovert->GetMean(),hyovert->GetRMS()); ttxt.DrawText(0.15,0.35, const_cast(txt)); p3->cd(); htminusy->DrawCopy(); mu0=mut-muy; sig0=sqrt(sigt*sigt+sigy*sigy); Ngauss->SetParameter(0, mu0); Ngauss->SetParameter(1, sig0); Ngauss->SetParameter(2, Nrndm * (xmax/2-xmin)/Nbins); Ngauss->DrawCopy("same"); p4->cd(); sprintf(txt,"T-Y Gauss: mu=%1.3g sig=%1.3g ==> Sd=%1.2g", mu0,sig0,fabs(mu0)/sig0); ttxt.DrawText(0.11,0.2, const_cast(txt)); sprintf(txt," ToyMC: mean=%1.3g RMS=%1.3g",htminusy->GetMean(),htminusy->GetRMS()); ttxt.DrawText(0.15,0.15, const_cast(txt)); // Save all objects in the ROOT file and close it. /* c1->Write(); hfile.Write(); hfile.Close(); */ // cout<< "weiter? "; cin >> int i; // delete c1; }