[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"
#include "TLinearFitter.h"

TObjArray arrayFit(3);
Int_t kMarkers[10]={25,24,20,21,22};
Int_t kColors[10]={1,2,4,3,6};
const char * names[10] = {"LTM","LThisto","LTMhistoPar","Gaus fit", "Gaus in range"};
const char * distNames[10] = {"Gauss","Gauss+flat","Gauss(m)+0.2*Gauss(m+5#sigma)","Gauss(m)+0.4*Gauss(m+5#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 meanLTMHistoPar(6);
  TVectorD sigmaLTMHistoPar(6);
  TVectorD meanLTMHistoGausLim(6);
  TVectorD sigmaLTMHistoGausLim(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.4) {                     //  gauss + second gaus 5 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];
      //
      // graph fit
      //
      Int_t binC=histo.GetXaxis()->FindBin(paramLTM[1]);
      Int_t bin0=TMath::Max(histo.GetXaxis()->FindBin(paramLTM[1]-5*paramLTM[2]),1);
      Int_t bin1=TMath::Min(histo.GetXaxis()->FindBin(paramLTM[1]+5*paramLTM[2]),histo.GetXaxis()->GetNbins()) ;
      if ( (bin0-bin1) <3){
	bin0=TMath::Max(binC-2,1);
	bin1=TMath::Min(binC+2, histo.GetXaxis()->GetNbins());
      }
      //
      /* test input:
	 TH1F histo("aaa","aaa",100,-10,10);for (Int_t i=0; i<10000; i++) histo.Fill(gRandom->Gaus(2,1));
      */
      
      TLinearFitter fitter(3,"hyp2");
      for (Int_t ibin=bin0; ibin<bin1; ibin++){
	Double_t x=histo.GetXaxis()->GetBinCenter(ibin);
	Double_t y=histo.GetBinContent(ibin);
	Double_t sy=histo.GetBinError(ibin);
	Double_t xxx[3]={x,x*x};
	if (y>3.*sy){
	  fitter.AddPoint(xxx,TMath::Log(y),sy/y);
	}
      }
      if (fitter.GetNpoints()>3){
	fitter.Eval();
	// f(x)=[0]+[1]*x+[2]*x*x;
	// f(x)=s*(x0-x0)^2+s
	Double_t parSigma=TMath::Sqrt(2*TMath::Abs(fitter.GetParameter(2)));
	Double_t parMean=-fitter.GetParameter(1)/(2.*fitter.GetParameter(2));  
	histo.Fit(&fg,"QN","QN", paramLTM[1]-5*paramLTM[2], paramLTM[1]+5*paramLTM[2]);
	meanLTMHistoPar[iltm]=parMean;
	sigmaLTMHistoPar[iltm]=parSigma;
	meanLTMHistoGausLim[iltm]=fg.GetParameter(1);
	sigmaLTMHistoGausLim[iltm]=fg.GetParameter(2);
      }else{
	meanLTMHistoPar[iltm]=paramLTM[1];
	sigmaLTMHistoPar[iltm]=paramLTM[2];
	meanLTMHistoGausLim[iltm]=meanG;
	sigmaLTMHistoGausLim[iltm]=rmsG;	
      }

    }
    (*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<<
      "meanLTMHistoPar.="<<&meanLTMHistoPar<<
      "meanLTMHistoGausLim.="<<&meanLTMHistoGausLim<<
      //
      "sigmaLTM.="<<&sigmaLTM<<
      "sigmaLTMHisto.="<<&sigmaLTMHisto<<
      "sigmaLTMHistoPar.="<<&sigmaLTMHistoPar<<      
      "sigmaLTMHistoGausLim.="<<&sigmaLTMHistoGausLim<<      
      "\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,0,0);
  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)*(meanLTMHistoPar.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();
    tree->Draw("sqrt(npoints-2)*(meanG-mean)/sigma:vecLTM.fElements>>hisLTMHist1(6,0.45,1.05,100,-10,10)",cutType+"npoints>50&&sigma>0.1","colzgoff");
    hisLTMTrunc[3]= (TH2*)tree->GetHistogram()->Clone();
    tree->Draw("sqrt(npoints-2)*(meanLTMHistoGausLim.fElements-mean)/sigma:vecLTM.fElements>>hisLTMHist1(6,0.45,1.05,100,-10,10)",cutType+"npoints>50&&sigma>0.1","colzgoff");
    hisLTMTrunc[4]= (TH2*)tree->GetHistogram()->Clone();
    TLegend * legend = new TLegend(0.5,0.7,0.88,0.88,distNames[itype]);
    legend->SetBorderSize(0);
    for (Int_t ihis=0; ihis<5; 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,1,kMarkers[ihis],kColors[ihis]);
      grSigma[ihis]->SetMaximum(5.0);
      grSigma[ihis]->SetMinimum(0.5);
      grSigma[ihis]->GetXaxis()->SetTitle("LTM fraction");
      grSigma[ihis]->GetYaxis()->SetTitle("#sigma_{meas}/#sigma_{Theor.}");
      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"
e5fec8fd	17	#include "TLinearFitter.h"
8ecd39c2	18
8ecd39c2	19	TObjArray arrayFit(3);
e5fec8fd	20	Int_t kMarkers[10]={25,24,20,21,22};
	21	Int_t kColors[10]={1,2,4,3,6};
	22	const char * names[10] = {"LTM","LThisto","LTMhistoPar","Gaus fit", "Gaus in range"};
	23	const char * distNames[10] = {"Gauss","Gauss+flat","Gauss(m)+0.2Gauss(m+5#sigma)","Gauss(m)+0.4Gauss(m+5#sigma)"};
8ecd39c2	24
	25
	26	void TestLTM(Int_t nevents=10000){
	27	//
	28	// Goal test and benchamerk numerical stability and precission of the LTM method
	29	// Binned data and not binned data:
	30	// Here we assume that binwidth<sigma
	31	// Distribution types examples used in test:
	32	//
	33	// 0 - gaussian
	34	// 1 - gauss+uniform background
	35	// 2 - gauss+second gaus 5 sigma away
	36	//
	37	Int_t npointsMax=1000;
	38	TTreeSRedirector * pcstream = new TTreeSRedirector("TStatToolkit_LTMtest.root","recreate");
	39	//
	40	TVectorD values(npointsMax);
	41	TVectorD vecLTM(6);
	42	TVectorD meanLTM(6);
	43	TVectorD sigmaLTM(6);
	44	TVectorD meanLTMHisto(6);
	45	TVectorD sigmaLTMHisto(6);
	46	TVectorD meanLTMHisto1(6);
	47	TVectorD sigmaLTMHisto1(6);
e5fec8fd	48	TVectorD meanLTMHistoPar(6);
	49	TVectorD sigmaLTMHistoPar(6);
	50	TVectorD meanLTMHistoGausLim(6);
	51	TVectorD sigmaLTMHistoGausLim(6);
8ecd39c2	52	TVectorD paramLTM(10);
	53	//
	54	for (Int_t iltm=0; iltm<6; iltm++) vecLTM[iltm]=0.99999*(0.5+iltm/10.);
	55	TF1 fg("fg","gaus");
	56	for (Int_t ievent = 0; ievent<nevents; ievent++){
	57	if (ievent%1000==0) printf("%d\n",ievent);
	58	Int_t distType=Int_t(gRandom->Rndm()*4);
	59	Int_t npoints= 50+(npointsMax-50)*gRandom->Rndm();
	60	Double_t mean = 0.5+(gRandom->Rndm()-0.5)*0.2;
	61	Double_t sigma = mean0.2(1+2*gRandom->Rndm());
	62	TH1F histo("histo","histo",100,-0.5,1.5);
	63	//
	64	for (Int_t ipoint=0; ipoint<npoints; ipoint++){
	65	Double_t value=0;
	66	if (distType==0) value=gRandom->Gaus(mean,sigma); // gauss
	67	if (distType==1) {
	68	if (gRandom->Rndm()>0.2) { // gauss + background
	69	value=gRandom->Gaus(mean,sigma);
	70	}else{
	71	value=mean+(gRandom->Rndm()-0.5)*2;
	72	}
	73	}
	74	if (distType==2) {
	75	if (gRandom->Rndm()>0.2) { // gauss + second gaus 5 sigma away
	76	value=gRandom->Gaus(mean,sigma);
	77	}else{
	78	value=gRandom->Gaus(mean+5*sigma,sigma);
	79	}
	80	}
	81	if (distType==3) {
e5fec8fd	82	if (gRandom->Rndm()>0.4) { // gauss + second gaus 5 sigma away
8ecd39c2	83	value=gRandom->Gaus(mean,sigma);
	84	}else{
	85	value=gRandom->Gaus(mean+5*sigma,sigma);
	86	}
	87	}
	88	values[ipoint]=value;
	89	histo.Fill(value);
	90	}
	91	//
	92	histo.Fit(&fg,"QN","QN");
	93	Double_t meanG = fg.GetParameter(1);
	94	Double_t rmsG = fg.GetParameter(2);
	95	Double_t meanA = TMath::Mean(npoints,values.GetMatrixArray());
	96	Double_t rmsA = TMath::Mean(npoints,values.GetMatrixArray());
	97	Double_t meanH = histo.GetMean();
	98	Double_t rmsH = histo.GetRMS();
	99	//
	100	for (Int_t iltm=0; iltm<6; iltm++){
	101	//
	102	Double_t meanV,sigmaV=0;
	103	TStatToolkit::EvaluateUni(npoints,values.GetMatrixArray(), meanV,sigmaV, vecLTM[iltm]*npoints);
	104	meanLTM[iltm]=meanV;
	105	sigmaLTM[iltm]=sigmaV;
	106	//
	107	TStatToolkit::LTM(&histo, &paramLTM, vecLTM[iltm]);
	108	meanLTMHisto1[iltm]=paramLTM[1];
	109	sigmaLTMHisto1[iltm]=paramLTM[2];
	110	//
	111	TStatToolkit::LTMHisto(&histo, paramLTM, vecLTM[iltm]);
	112	meanLTMHisto[iltm]=paramLTM[1];
	113	sigmaLTMHisto[iltm]=paramLTM[2];
e5fec8fd	114	//
	115	// graph fit
	116	//
	117	Int_t binC=histo.GetXaxis()->FindBin(paramLTM[1]);
	118	Int_t bin0=TMath::Max(histo.GetXaxis()->FindBin(paramLTM[1]-5*paramLTM[2]),1);
	119	Int_t bin1=TMath::Min(histo.GetXaxis()->FindBin(paramLTM[1]+5*paramLTM[2]),histo.GetXaxis()->GetNbins()) ;
	120	if ( (bin0-bin1) <3){
	121	bin0=TMath::Max(binC-2,1);
	122	bin1=TMath::Min(binC+2, histo.GetXaxis()->GetNbins());
	123	}
	124	//
	125	/* test input:
	126	TH1F histo("aaa","aaa",100,-10,10);for (Int_t i=0; i<10000; i++) histo.Fill(gRandom->Gaus(2,1));
	127	*/
	128
	129	TLinearFitter fitter(3,"hyp2");
	130	for (Int_t ibin=bin0; ibin<bin1; ibin++){
	131	Double_t x=histo.GetXaxis()->GetBinCenter(ibin);
	132	Double_t y=histo.GetBinContent(ibin);
	133	Double_t sy=histo.GetBinError(ibin);
	134	Double_t xxx[3]={x,x*x};
	135	if (y>3.*sy){
	136	fitter.AddPoint(xxx,TMath::Log(y),sy/y);
	137	}
	138	}
	139	if (fitter.GetNpoints()>3){
	140	fitter.Eval();
	141	// f(x)=[0]+[1]x+[2]x*x;
	142	// f(x)=s*(x0-x0)^2+s
	143	Double_t parSigma=TMath::Sqrt(2*TMath::Abs(fitter.GetParameter(2)));
	144	Double_t parMean=-fitter.GetParameter(1)/(2.*fitter.GetParameter(2));
	145	histo.Fit(&fg,"QN","QN", paramLTM[1]-5paramLTM[2], paramLTM[1]+5paramLTM[2]);
	146	meanLTMHistoPar[iltm]=parMean;
	147	sigmaLTMHistoPar[iltm]=parSigma;
	148	meanLTMHistoGausLim[iltm]=fg.GetParameter(1);
	149	sigmaLTMHistoGausLim[iltm]=fg.GetParameter(2);
	150	}else{
	151	meanLTMHistoPar[iltm]=paramLTM[1];
	152	sigmaLTMHistoPar[iltm]=paramLTM[2];
	153	meanLTMHistoGausLim[iltm]=meanG;
	154	sigmaLTMHistoGausLim[iltm]=rmsG;
	155	}
	156
8ecd39c2	157	}
	158	(*pcstream)<<"ltm"<<
	159	//Input
	160	"npoints="<<npoints<<
	161	"distType="<<distType<<
	162	"mean="<<mean<<
	163	"sigma="<<sigma<<
	164	// "Standart" statistic output
	165	"meanA="<<meanA<<
	166	"rmsA="<<rmsA<<
	167	"meanH="<<meanH<<
	168	"rmsH="<<rmsH<<
	169	"meanG="<<meanG<<
	170	"rmsG="<<rmsG<<
	171	// "LTM" output
	172	"vecLTM.="<<&vecLTM<<
	173	"meanLTM.="<<&meanLTM<<
	174	"meanLTMHisto.="<<&meanLTMHisto<<
	175	"meanLTMHisto1.="<<&meanLTMHisto1<<
e5fec8fd	176	"meanLTMHistoPar.="<<&meanLTMHistoPar<<
e5fec8fd	177	"meanLTMHistoGausLim.="<<&meanLTMHistoGausLim<<
8ecd39c2	178	//
	179	"sigmaLTM.="<<&sigmaLTM<<
	180	"sigmaLTMHisto.="<<&sigmaLTMHisto<<
e5fec8fd	181	"sigmaLTMHistoPar.="<<&sigmaLTMHistoPar<<
e5fec8fd	182	"sigmaLTMHistoGausLim.="<<&sigmaLTMHistoGausLim<<
8ecd39c2	183	"\n";
	184	}
	185	delete pcstream;
	186	//
	187	//
	188	//
	189	TFile *fltm = TFile::Open("TStatToolkit_LTMtest.root","update");
	190	TTree * tree = (TTree*)fltm->Get("ltm");
	191	tree->SetMarkerSize(0.5);
	192	tree->SetMarkerStyle(25);
	193	//
	194	// 1. Get numerical error estimate of the LTM method for gaussian distribution
	195	//
	196	TH2 * hisLTMTrunc[10]={0};
	197	TH1 * hisResol[10]={0};
	198	TGraphErrors * grSigma[10]={0};
	199	TCanvas *canvasLTM= new TCanvas("canvasLTM","canvasLTM",800,700);
e5fec8fd	200	canvasLTM->Divide(2,2,0,0);
8ecd39c2	201	for (Int_t itype=0; itype<4; itype++){
	202	canvasLTM->cd(itype+1);
	203	TCut cutType=TString::Format("distType==0%d",itype).Data();
	204	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");
	205	hisLTMTrunc[0]= (TH2*)(tree->GetHistogram()->Clone());
	206	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");
	207	hisLTMTrunc[1]= (TH2*)tree->GetHistogram()->Clone();
e5fec8fd	208	tree->Draw("sqrt(npoints-2)*(meanLTMHistoPar.fElements-mean)/sigma:vecLTM.fElements>>hisLTMHist1(6,0.45,1.05,100,-10,10)",cutType+"npoints>50&&sigma>0.1","colzgoff");
8ecd39c2	209	hisLTMTrunc[2]= (TH2*)tree->GetHistogram()->Clone();
e5fec8fd	210	tree->Draw("sqrt(npoints-2)*(meanG-mean)/sigma:vecLTM.fElements>>hisLTMHist1(6,0.45,1.05,100,-10,10)",cutType+"npoints>50&&sigma>0.1","colzgoff");
	211	hisLTMTrunc[3]= (TH2*)tree->GetHistogram()->Clone();
	212	tree->Draw("sqrt(npoints-2)*(meanLTMHistoGausLim.fElements-mean)/sigma:vecLTM.fElements>>hisLTMHist1(6,0.45,1.05,100,-10,10)",cutType+"npoints>50&&sigma>0.1","colzgoff");
	213	hisLTMTrunc[4]= (TH2*)tree->GetHistogram()->Clone();
	214	TLegend * legend = new TLegend(0.5,0.7,0.88,0.88,distNames[itype]);
8ecd39c2	215	legend->SetBorderSize(0);
e5fec8fd	216	for (Int_t ihis=0; ihis<5; ihis++){
8ecd39c2	217	// MakeStat1D(TH2 * his, Int_t deltaBin, Double_t fraction, Int_t returnType, Int_t markerStyle, Int_t markerColor);
e5fec8fd	218	grSigma[ihis]=TStatToolkit::MakeStat1D( hisLTMTrunc[ihis],0,0.99,1,kMarkers[ihis],kColors[ihis]);
	219	grSigma[ihis]->SetMaximum(5.0);
	220	grSigma[ihis]->SetMinimum(0.5);
8ecd39c2	221	grSigma[ihis]->GetXaxis()->SetTitle("LTM fraction");
e5fec8fd	222	grSigma[ihis]->GetYaxis()->SetTitle("#sigma_{meas}/#sigma_{Theor.}");
8ecd39c2	223	if (ihis==0)grSigma[ihis]->Draw("alp");
	224	if (ihis>0)grSigma[ihis]->Draw("lp");
	225	legend->AddEntry(grSigma[ihis],names[ihis],"p");
	226	}
	227	legend->Draw();
	228	}
	229	canvasLTM->SaveAs("robustStatLTM_Performance.pdf");
	230
	231
	232	}