// EF evaluates the systematic error in BR for sigma t0

#include "TH1.h"
#include "TF1.h"
#include "TRandom3.h"

#include <fstream>
#include <iostream>

#include <stdio.h>
#include <stdlib.h>

// Output text files names
ifstream in("his.dat");
ofstream out("out.dat");

Double_t t_mich(Double_t *t, Double_t *par) // single Michel time pdf
 { 
   double yield = 0.;
   //par[1]=0;
   double tau_pi = 26.03;
   double tau_mu = 2197.019;
   double pg_start = -100.; 
   double pg_close = 200;

   if ( t[0]>par[1] ) { yield = par[0]*(exp(-(t[0]+par[1])/tau_mu) - exp(-(t[0]+par[1])/tau_pi)); return yield; }
   if ( t[0] <=par[1] ) return par[0]*yield;
}

void pgate(Int_t nrun=1)
 {

  TRandom3 *r = new TRandom3();

  //TRandom3 r;
  r->SetSeed(0);

   double frac[20], delta[20];

  TFile *file  = new TFile("histos.root","RECREATE");

  //TF1 *t_mich = new TF1("t_mich","exp(-x/2197.019)-exp(-x/26.03)",0,200);
  //  TF1 *t_mich = new TF1("t_mich",t_mich,-100,200,2);
  TF1 *t_mich = new TF1("t_mich",t_mich,60,140,2);

  TH1F *h900010 = new TH1F("h900010","simulated michel decay time",160, 60,140);
  TH1F *h900011 = new TH1F("h900011","simulated michel decay time",160, 60,140);
  TH1F *h900012 = new TH1F("h900012","simulated michel decay time",160, 60,140);

  Double_t par[2]={1.,0};
 t_mich->SetParameters(par);

  t_mich->Draw();

  /*
  for (int i=-10;i<10;i++) {
    float x=200+ 0.02*i;
    j=10+i;
    frac[j]=t_mich->Integral(0,x)/t_mich->Integral(0,50000);
    printf(" time: %f fraction: %f\n",x,t_mich->Integral(0,x)/t_mich->Integral(0,100000));
  }

     for (int i=0;i<19;i++) {
       delta[i]=100*(frac[i+1]-frac[i])/frac[i];
       printf(" delta pecrent %f\n",delta[i]);
     }
  */

   // create random histogram

  //TH1D *h1 = new TH1D("h1","h1",300,-100.,200.);

  //h1->FillRandom("t_mich",1E8);

  //h1->Fit(t_mich,"","",60,200);  
  //h1->Fit("pol5","","",150,200); 

  //h900012->FillRandom("t_mich",600);
  //r->SetSeed(0);
  h900010->FillRandom("t_mich",606300);

    printf (" f200= %f\n",t_mich->Eval(200.));
  // printf (" f200= %f\n",pol5->Eval(200.));

   h900010->SetMarkerStyle(20);
   h900010->SetMarkerColor(kRed);
   h900010->SetMarkerSize(0.5);
   h900011->SetMarkerStyle(20);
   h900011->SetMarkerColor(kRed);
   h900011->SetMarkerSize(0.5);

   h900011->Add(h900010,h900010,0.5,0.5);  

   c1 = new TCanvas("c1"," ",10,10,700,1000);
   c1->Clear();
   c1->Divide(1,2);
   c1->cd(1);
   h900010->Draw("e1");
   c1->cd(2);
   h900011->Smooth(10);
   h900011->Draw();
   c1->Update();
   c1->Print("michel_r1.eps");

   file->Write();

}