STAT/Macros/TStatToolkitTest.C

   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 #include "TLinearFitter.h"
  18
  19 TObjArray arrayFit(3);
  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.2*Gauss(m+5#sigma)","Gauss(m)+0.4*Gauss(m+5#sigma)"};
  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);
  48   TVectorD meanLTMHistoPar(6);
  49   TVectorD sigmaLTMHistoPar(6);
  50   TVectorD meanLTMHistoGausLim(6);
  51   TVectorD sigmaLTMHistoGausLim(6);
  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 = mean*0.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) {
  82         if (gRandom->Rndm()>0.4) {                     //  gauss + second gaus 5 sigma away
  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];
 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]-5*paramLTM[2], paramLTM[1]+5*paramLTM[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
 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<<
 176       "meanLTMHistoPar.="<<&meanLTMHistoPar<<
 177       "meanLTMHistoGausLim.="<<&meanLTMHistoGausLim<<
 178       //
 179       "sigmaLTM.="<<&sigmaLTM<<
 180       "sigmaLTMHisto.="<<&sigmaLTMHisto<<
 181       "sigmaLTMHistoPar.="<<&sigmaLTMHistoPar<<
 182       "sigmaLTMHistoGausLim.="<<&sigmaLTMHistoGausLim<<
 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);
 200   canvasLTM->Divide(2,2,0,0);
 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();
 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");
 209     hisLTMTrunc[2]= (TH2*)tree->GetHistogram()->Clone();
 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]);
 215     legend->SetBorderSize(0);
 216     for (Int_t ihis=0; ihis<5; ihis++){
 217       // MakeStat1D(TH2 * his, Int_t deltaBin, Double_t fraction, Int_t returnType, Int_t markerStyle, Int_t markerColor);
 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);
 221       grSigma[ihis]->GetXaxis()->SetTitle("LTM fraction");
 222       grSigma[ihis]->GetYaxis()->SetTitle("#sigma_{meas}/#sigma_{Theor.}");
 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 }