[u/mrichter/AliRoot.git] / STAT / Macros / TStatToolkitTest.C

/*
  .L $ALICE_ROOT/STAT/Macros/TStatToolkitTest.C+

*/

#include "TH1.h"
#include "TF1.h"
#include "TMath.h"
#include "TRandom.h"
#include "TVectorD.h"
#include "TStatToolkit.h"
#include "TTreeStream.h"
#include "TLegend.h"
#include "TGraphErrors.h"
#include "TCut.h"
#include "TCanvas.h"

TObjArray arrayFit(3);
Int_t kMarkers[10]={25,24,20,21};
Int_t kColors[10]={1,2,4,3};
const char * names[10] = {"LTM","LThisto","LTMhisto1"};
const char * distNames[10] = {"Gauss","Gauss+flat","Gauss(m)+0.2*Gauss(m+5#sigma)","Gauss(m)+0.3*Gauss(m+3#sigma)"};


void TestLTM(Int_t nevents=10000){
  //
  // Goal test and benchamerk numerical stability and precission of the LTM method
  // Binned data and not binned data: 
  // Here we assume that binwidth<sigma
  // Distribution types examples used in test:
  //
  //   0 - gaussian
  //   1 - gauss+uniform background
  //   2 - gauss+second gaus 5 sigma away 
  //   
  Int_t npointsMax=1000;
  TTreeSRedirector * pcstream = new TTreeSRedirector("TStatToolkit_LTMtest.root","recreate");
  //
  TVectorD values(npointsMax);
  TVectorD vecLTM(6);
  TVectorD meanLTM(6);
  TVectorD sigmaLTM(6);
  TVectorD meanLTMHisto(6);
  TVectorD sigmaLTMHisto(6);
  TVectorD meanLTMHisto1(6);
  TVectorD sigmaLTMHisto1(6);
  TVectorD paramLTM(10);
  //
  for (Int_t iltm=0; iltm<6; iltm++) vecLTM[iltm]=0.99999*(0.5+iltm/10.);
  TF1 fg("fg","gaus");
  for (Int_t ievent = 0; ievent<nevents; ievent++){
    if (ievent%1000==0) printf("%d\n",ievent);
    Int_t distType=Int_t(gRandom->Rndm()*4);
    Int_t npoints= 50+(npointsMax-50)*gRandom->Rndm();
    Double_t mean  = 0.5+(gRandom->Rndm()-0.5)*0.2;
    Double_t sigma = mean*0.2*(1+2*gRandom->Rndm());
    TH1F histo("histo","histo",100,-0.5,1.5);
    //
    for (Int_t ipoint=0; ipoint<npoints; ipoint++){
      Double_t value=0;
      if (distType==0) value=gRandom->Gaus(mean,sigma); // gauss
      if (distType==1) {
	if (gRandom->Rndm()>0.2) {                     //  gauss + background
	  value=gRandom->Gaus(mean,sigma);
	}else{
	  value=mean+(gRandom->Rndm()-0.5)*2;
	}
      }
      if (distType==2) {
	if (gRandom->Rndm()>0.2) {                     //  gauss + second gaus 5 sigma away
	  value=gRandom->Gaus(mean,sigma);
	}else{
	  value=gRandom->Gaus(mean+5*sigma,sigma);
	}
      }
      if (distType==3) {
	if (gRandom->Rndm()>0.3) {                     //  gauss + second gaus 4 sigma away
	  value=gRandom->Gaus(mean,sigma);
	}else{
	  value=gRandom->Gaus(mean+5*sigma,sigma);
	}
      }
      values[ipoint]=value;
      histo.Fill(value);
    }
    //
    histo.Fit(&fg,"QN","QN");
    Double_t meanG = fg.GetParameter(1);
    Double_t rmsG = fg.GetParameter(2);
    Double_t meanA  = TMath::Mean(npoints,values.GetMatrixArray());
    Double_t rmsA   = TMath::Mean(npoints,values.GetMatrixArray());
    Double_t meanH  = histo.GetMean();
    Double_t rmsH  = histo.GetRMS();
    //
    for (Int_t iltm=0; iltm<6; iltm++){
      //    
      Double_t meanV,sigmaV=0;
      TStatToolkit::EvaluateUni(npoints,values.GetMatrixArray(), meanV,sigmaV, vecLTM[iltm]*npoints);
      meanLTM[iltm]=meanV;
      sigmaLTM[iltm]=sigmaV;
      //
      TStatToolkit::LTM(&histo, &paramLTM,  vecLTM[iltm]);
      meanLTMHisto1[iltm]=paramLTM[1];
      sigmaLTMHisto1[iltm]=paramLTM[2];
      //
      TStatToolkit::LTMHisto(&histo, paramLTM,  vecLTM[iltm]);
      meanLTMHisto[iltm]=paramLTM[1];
      sigmaLTMHisto[iltm]=paramLTM[2];
    }
    (*pcstream)<<"ltm"<<
      //Input
      "npoints="<<npoints<<
      "distType="<<distType<<
      "mean="<<mean<<
      "sigma="<<sigma<<
      // "Standart" statistic output
      "meanA="<<meanA<<
      "rmsA="<<rmsA<<
      "meanH="<<meanH<<
      "rmsH="<<rmsH<<
      "meanG="<<meanG<<
      "rmsG="<<rmsG<<
      // "LTM" output
      "vecLTM.="<<&vecLTM<<
      "meanLTM.="<<&meanLTM<<
      "meanLTMHisto.="<<&meanLTMHisto<<
      "meanLTMHisto1.="<<&meanLTMHisto1<<
      //
      "sigmaLTM.="<<&sigmaLTM<<
      "sigmaLTMHisto.="<<&sigmaLTMHisto<<
      "sigmaLTMHisto1.="<<&sigmaLTMHisto1<<      
      "\n";
  }
  delete pcstream;
  //
  //
  //
  TFile *fltm = TFile::Open("TStatToolkit_LTMtest.root","update");
  TTree * tree = (TTree*)fltm->Get("ltm");
  tree->SetMarkerSize(0.5);
  tree->SetMarkerStyle(25);
  //
  // 1. Get numerical error estimate of the LTM method for gaussian distribution  
  //
  TH2 * hisLTMTrunc[10]={0};  
  TH1 * hisResol[10]={0};
  TGraphErrors * grSigma[10]={0};
  TCanvas *canvasLTM= new TCanvas("canvasLTM","canvasLTM",800,700);
  canvasLTM->Divide(2,2);
  for (Int_t itype=0; itype<4; itype++){
    canvasLTM->cd(itype+1);
    TCut cutType=TString::Format("distType==0%d",itype).Data();
    tree->Draw("sqrt(npoints-2)*(meanLTM.fElements-mean)/sigma:vecLTM.fElements>>hisLTM(6,0.45,1.05,100,-10,10)",cutType+"npoints>50&&sigma>0.1","colzgof");
    hisLTMTrunc[0]= (TH2*)(tree->GetHistogram()->Clone());
    tree->Draw("sqrt(npoints-2)*(meanLTMHisto.fElements-mean)/sigma:vecLTM.fElements>>hisLTMHisto(6,0.45,1.05,100,-10,10)",cutType+"npoints>50&&sigma>0.1","colzgoff");
    hisLTMTrunc[1]= (TH2*)tree->GetHistogram()->Clone();
    tree->Draw("sqrt(npoints-2)*(meanLTMHisto1.fElements-mean)/sigma:vecLTM.fElements>>hisLTMHist1(6,0.45,1.05,100,-10,10)",cutType+"npoints>50&&sigma>0.1","colzgoff");
    hisLTMTrunc[2]= (TH2*)tree->GetHistogram()->Clone();
    TLegend * legend = new TLegend(0.5,0.7,0.9,0.9,distNames[itype]);
    legend->SetBorderSize(0);
    for (Int_t ihis=0; ihis<3; ihis++){
      // MakeStat1D(TH2 * his, Int_t deltaBin, Double_t fraction, Int_t returnType, Int_t markerStyle, Int_t markerColor);
      grSigma[ihis]=TStatToolkit::MakeStat1D( hisLTMTrunc[ihis],0,0.99,5,kMarkers[ihis],kColors[ihis]);
      grSigma[ihis]->GetXaxis()->SetTitle("LTM fraction");
      grSigma[ihis]->GetYaxis()->SetTitle("#sigma_{meas}/sigma_{gauss}");
      if (ihis==0)grSigma[ihis]->Draw("alp");
      if (ihis>0)grSigma[ihis]->Draw("lp");
      legend->AddEntry(grSigma[ihis],names[ihis],"p");
    }
    legend->Draw();
  }
  canvasLTM->SaveAs("robustStatLTM_Performance.pdf");


}
Commit	Line	Data
8ecd39c2	1	/*
	2	.L $ALICE_ROOT/STAT/Macros/TStatToolkitTest.C+
	3
	4	*/
	5
	6	#include "TH1.h"
	7	#include "TF1.h"
	8	#include "TMath.h"
	9	#include "TRandom.h"
	10	#include "TVectorD.h"
	11	#include "TStatToolkit.h"
	12	#include "TTreeStream.h"
	13	#include "TLegend.h"
	14	#include "TGraphErrors.h"
	15	#include "TCut.h"
	16	#include "TCanvas.h"
	17
	18	TObjArray arrayFit(3);
	19	Int_t kMarkers[10]={25,24,20,21};
	20	Int_t kColors[10]={1,2,4,3};
	21	const char * names[10] = {"LTM","LThisto","LTMhisto1"};
	22	const char * distNames[10] = {"Gauss","Gauss+flat","Gauss(m)+0.2Gauss(m+5#sigma)","Gauss(m)+0.3Gauss(m+3#sigma)"};
	23
	24
	25	void TestLTM(Int_t nevents=10000){
	26	//
	27	// Goal test and benchamerk numerical stability and precission of the LTM method
	28	// Binned data and not binned data:
	29	// Here we assume that binwidth<sigma
	30	// Distribution types examples used in test:
	31	//
	32	// 0 - gaussian
	33	// 1 - gauss+uniform background
	34	// 2 - gauss+second gaus 5 sigma away
	35	//
	36	Int_t npointsMax=1000;
	37	TTreeSRedirector * pcstream = new TTreeSRedirector("TStatToolkit_LTMtest.root","recreate");
	38	//
	39	TVectorD values(npointsMax);
	40	TVectorD vecLTM(6);
	41	TVectorD meanLTM(6);
	42	TVectorD sigmaLTM(6);
	43	TVectorD meanLTMHisto(6);
	44	TVectorD sigmaLTMHisto(6);
	45	TVectorD meanLTMHisto1(6);
	46	TVectorD sigmaLTMHisto1(6);
	47	TVectorD paramLTM(10);
	48	//
	49	for (Int_t iltm=0; iltm<6; iltm++) vecLTM[iltm]=0.99999*(0.5+iltm/10.);
	50	TF1 fg("fg","gaus");
	51	for (Int_t ievent = 0; ievent<nevents; ievent++){
	52	if (ievent%1000==0) printf("%d\n",ievent);
	53	Int_t distType=Int_t(gRandom->Rndm()*4);
	54	Int_t npoints= 50+(npointsMax-50)*gRandom->Rndm();
	55	Double_t mean = 0.5+(gRandom->Rndm()-0.5)*0.2;
	56	Double_t sigma = mean0.2(1+2*gRandom->Rndm());
	57	TH1F histo("histo","histo",100,-0.5,1.5);
	58	//
	59	for (Int_t ipoint=0; ipoint<npoints; ipoint++){
	60	Double_t value=0;
	61	if (distType==0) value=gRandom->Gaus(mean,sigma); // gauss
	62	if (distType==1) {
	63	if (gRandom->Rndm()>0.2) { // gauss + background
	64	value=gRandom->Gaus(mean,sigma);
65	}else{
66	value=mean+(gRandom->Rndm()-0.5)*2;
67	}
68	}
69	if (distType==2) {
70	if (gRandom->Rndm()>0.2) { // gauss + second gaus 5 sigma away
71	value=gRandom->Gaus(mean,sigma);
72	}else{
73	value=gRandom->Gaus(mean+5*sigma,sigma);
74	}
75	}
76	if (distType==3) {
77	if (gRandom->Rndm()>0.3) { // gauss + second gaus 4 sigma away
78	value=gRandom->Gaus(mean,sigma);
79	}else{
80	value=gRandom->Gaus(mean+5*sigma,sigma);
81	}
82	}
83	values[ipoint]=value;
84	histo.Fill(value);
85	}
86	//
87	histo.Fit(&fg,"QN","QN");
88	Double_t meanG = fg.GetParameter(1);
89	Double_t rmsG = fg.GetParameter(2);
90	Double_t meanA = TMath::Mean(npoints,values.GetMatrixArray());
91	Double_t rmsA = TMath::Mean(npoints,values.GetMatrixArray());
92	Double_t meanH = histo.GetMean();
93	Double_t rmsH = histo.GetRMS();
94	//
95	for (Int_t iltm=0; iltm<6; iltm++){
96	//
97	Double_t meanV,sigmaV=0;
98	TStatToolkit::EvaluateUni(npoints,values.GetMatrixArray(), meanV,sigmaV, vecLTM[iltm]*npoints);
99	meanLTM[iltm]=meanV;
100	sigmaLTM[iltm]=sigmaV;
101	//
102	TStatToolkit::LTM(&histo, &paramLTM, vecLTM[iltm]);
103	meanLTMHisto1[iltm]=paramLTM[1];
104	sigmaLTMHisto1[iltm]=paramLTM[2];
105	//
106	TStatToolkit::LTMHisto(&histo, paramLTM, vecLTM[iltm]);
107	meanLTMHisto[iltm]=paramLTM[1];
108	sigmaLTMHisto[iltm]=paramLTM[2];
109	}
110	(*pcstream)<<"ltm"<<
111	//Input
112	"npoints="<<npoints<<
113	"distType="<<distType<<
114	"mean="<<mean<<
115	"sigma="<<sigma<<
116	// "Standart" statistic output
117	"meanA="<<meanA<<
118	"rmsA="<<rmsA<<
119	"meanH="<<meanH<<
120	"rmsH="<<rmsH<<
121	"meanG="<<meanG<<
122	"rmsG="<<rmsG<<
123	// "LTM" output
124	"vecLTM.="<<&vecLTM<<
125	"meanLTM.="<<&meanLTM<<
126	"meanLTMHisto.="<<&meanLTMHisto<<
127	"meanLTMHisto1.="<<&meanLTMHisto1<<
128	//
129	"sigmaLTM.="<<&sigmaLTM<<
130	"sigmaLTMHisto.="<<&sigmaLTMHisto<<
131	"sigmaLTMHisto1.="<<&sigmaLTMHisto1<<
132	"\n";
133	}
134	delete pcstream;
135	//
136	//
137	//
138	TFile *fltm = TFile::Open("TStatToolkit_LTMtest.root","update");
139	TTree * tree = (TTree*)fltm->Get("ltm");
140	tree->SetMarkerSize(0.5);
141	tree->SetMarkerStyle(25);
142	//
143	// 1. Get numerical error estimate of the LTM method for gaussian distribution
144	//
145	TH2 * hisLTMTrunc[10]={0};
146	TH1 * hisResol[10]={0};
147	TGraphErrors * grSigma[10]={0};
148	TCanvas *canvasLTM= new TCanvas("canvasLTM","canvasLTM",800,700);
149	canvasLTM->Divide(2,2);
150	for (Int_t itype=0; itype<4; itype++){
151	canvasLTM->cd(itype+1);
152	TCut cutType=TString::Format("distType==0%d",itype).Data();
153	tree->Draw("sqrt(npoints-2)*(meanLTM.fElements-mean)/sigma:vecLTM.fElements>>hisLTM(6,0.45,1.05,100,-10,10)",cutType+"npoints>50&&sigma>0.1","colzgof");
154	hisLTMTrunc[0]= (TH2*)(tree->GetHistogram()->Clone());
155	tree->Draw("sqrt(npoints-2)*(meanLTMHisto.fElements-mean)/sigma:vecLTM.fElements>>hisLTMHisto(6,0.45,1.05,100,-10,10)",cutType+"npoints>50&&sigma>0.1","colzgoff");
156	hisLTMTrunc[1]= (TH2*)tree->GetHistogram()->Clone();
157	tree->Draw("sqrt(npoints-2)*(meanLTMHisto1.fElements-mean)/sigma:vecLTM.fElements>>hisLTMHist1(6,0.45,1.05,100,-10,10)",cutType+"npoints>50&&sigma>0.1","colzgoff");
158	hisLTMTrunc[2]= (TH2*)tree->GetHistogram()->Clone();
159	TLegend * legend = new TLegend(0.5,0.7,0.9,0.9,distNames[itype]);
160	legend->SetBorderSize(0);
161	for (Int_t ihis=0; ihis<3; ihis++){
162	// MakeStat1D(TH2 * his, Int_t deltaBin, Double_t fraction, Int_t returnType, Int_t markerStyle, Int_t markerColor);
163	grSigma[ihis]=TStatToolkit::MakeStat1D( hisLTMTrunc[ihis],0,0.99,5,kMarkers[ihis],kColors[ihis]);
164	grSigma[ihis]->GetXaxis()->SetTitle("LTM fraction");
165	grSigma[ihis]->GetYaxis()->SetTitle("#sigma_{meas}/sigma_{gauss}");
166	if (ihis==0)grSigma[ihis]->Draw("alp");
167	if (ihis>0)grSigma[ihis]->Draw("lp");
168	legend->AddEntry(grSigma[ihis],names[ihis],"p");
169	}
170	legend->Draw();
171	}
172	canvasLTM->SaveAs("robustStatLTM_Performance.pdf");
173
174
175	}