void test_MeanError() {

   // test mean error in case of weighted histograms

   int n1 = 500;
   int n2 = 400; 
   auto h1 = new TH1D("h1","h1",100,-3,3); 
   auto h2 = new TH1D("h2","h2",100,-3,3);
   auto hdiff =  new TH1D("hdiff","h1-h2",100,-3,3);
   h1->Sumw2();
   h2->Sumw2();

   auto hmean = new TH1D("hmean","mean ",50,1,0);
   auto hpull = new TH1D("hpull","pull ",50,-3,3); 
   auto hpull2 = new TH1D("hpull2","pull2 ",50,-3,3); 

   int ntoys = 10000;

   for (int ievt = 0; ievt < ntoys; ++ievt) { 

      h1->Reset();
      h2->Reset();
      hdiff->Reset();
      
      h1->FillRandom("gaus",n1); 
      h2->FillRandom("gaus",n2); 

      hdiff->Add(h1,h2, 1, -1);

      hmean->Fill(hdiff->GetMean());
      hpull->Fill(hdiff->GetMean()/hdiff->GetMeanError() );

      // use formula from Cocrhan (1977) and reference in this paper
      // https://www.sciencedirect.com/science/article/pii/135223109400210C

      double s[10];
      hdiff->GetStats(s); 
      
      double n = hdiff->GetNbinsX();
      double xWbar = hdiff->GetMean();
      double sumw = s[0]; 
      double wbar = sumw/double(n);   // s[0] is sumw

      // formula 
      // out = n/((n-1)*sum(w)^2)*(sum((w*x-wbar*xWbar)^2)-2*xWbar*sum((w-wbar)*(w*x-wbar*xWbar))+xWbar^2*sum((w-wbar)^2))
      double sum1 = 0;
      double sum2 = 0;
      double sum3 = 0; 

      for (int i =0; i < n; ++i ) {
         double w = hdiff->GetBinContent(i+1);
         double x = hdiff->GetBinCenter(i+1); 
         double tmp =   w*x - wbar*xWbar;
         sum1 += tmp*tmp;

         sum2 += ( w - wbar)*(w*x-wbar*xWbar);

         sum3 += ( w - wbar)*(w-wbar);
      }
         
      double result  = double(n)/(double(n-1)*sumw*sumw ) * ( sum1 - 2.*xWbar * sum2 + xWbar*xWbar * sum3 );

      hpull2->Fill( xWbar / sqrt(result) ); 

      
   } 

   auto c1 = new TCanvas();
   c1->Divide(1,3);
   c1->cd(1);
   hmean->Draw();
   c1->cd(2);
   hpull->Draw();
   hpull->Fit("gaus","L"); 
   c1->cd(3);
   hpull2->Draw(); 
   hpull2->Fit("gaus","L"); 


   // compute chi2 after fit
   std::cout << "chi2 for pull 1 " <<   hpull->Chisquare(hpull->GetFunction("gaus"), "L") << std::endl; 
   std::cout << "chi2 for pull 2 " <<   hpull2->Chisquare(hpull2->GetFunction("gaus"), "L") << std::endl; 
}