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
  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.2*Gauss(m+5#sigma)","Gauss(m)+0.3*Gauss(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 = mean*0.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 }