using namespace TMath;

const double pi = Pi();

const double xmin = -0.01;
const double xmax =  0.01;

enum kPars  { kfracc, kfrac25, kfrac50, 
               kmuc, ksc, 
               kmu25, ks25, 
               kmu50, ks50, 
               kmuO, ksO,
               knorm 
             };

double g( double *x, double *par ) {

  double xx = x[0];
  double mu = par[0];
  double s  = par[1];

  double t  = xx-mu;
         t  /= s;
         t  *= t;
         t  /= 2.;

  double norm = Sqrt(2.*pi)*s;

  double exp  = Exp( - t );

  return exp/norm;

}

double func( double *x, double *par ) {

  double xx = x[0];

  double fracc   = par[kfracc];
  double frac25  = par[kfrac25];
  double frac50  = par[kfrac50];


  double muc     = par[kmuc];
  double sc     = par[ksc];

  double mu25    = par[kmu25];
  double s25    = par[ks25];

  double mu50    = par[kmu50];
  double s50    = par[ks50];

  double muO     = par[kmuO];
  double sO     = par[ksO]; 

  double normaliz = par[knorm];

  double parc[2]  = { muc, sc };

  double par25[2] = { mu25, s25 };

  double par50[2] = { mu50, s50 };

  double parO[2] =  { muO, sO };

  double value  =  fracc * (g(x,parc));

         value  += (1.-fracc)*(frac25)*g(x,par25);
         par25[0] = -mu25;
         value  += (1.-fracc)*(frac25)*g(x,par25);

         value  += (1.-fracc)*(1.-frac25)*(frac50)*g(x,par50);
         par50[0] = -mu50;
         value  += (1.-fracc)*(1.-frac25)*(frac50)*g(x,par50);

         value  += (1.-fracc)*(1.-frac25)*(1.-frac50)*g(x,parO);

  return normaliz*value;
}

void bad_res(Int_t runNumber) {

  TString fileName(Form("rootuple/Sbt_Run%d_resout.root",runNumber));
  TFile file(fileName.Data(),"READ");
  file.cd();
  TH1D* uRes = (TH1D*)file.Get("ResU_bad_Det4");
  TH1D* vRes = (TH1D*)file.Get("ResV_bad_Det4");
  

  TF1 f("f",func,xmin,xmax,12);

  f.SetParName(kfracc,"fracc");
  f.SetParName(kfrac25,"frac25");
  f.SetParName(kfrac50,"frac50");

  f.SetParName(kmuc,"muc");
  f.SetParName(ksc,"sc");

  f.SetParName(kmu25,"mu25");
  f.SetParName(ks25,"s25");

  f.SetParName(kmu50,"mu50");
  f.SetParName(ks50,"s50");

  f.SetParName(kmuO,"muO");
  f.SetParName(ksO,"sO");

  f.SetParName(knorm,"norm");


  f.SetParLimits(kfracc,.0,1.);  
  f.SetParLimits(kfrac25,.0,1.);  
  f.SetParLimits(kfrac50,.0,1.);  

  f.SetParameter(kfracc,.7);
  f.SetParameter(kfrac25,.5);
  f.SetParameter(kfrac50,.5);

  
  f.SetParLimits(kmuc,-1.e-3,1.e-3);
  f.SetParLimits(ksc,+1.e-3,3.e-3);
  
  f.SetParameter(kmuc,0.);
  f.SetParameter(ksc,2.e-3);

  
  f.SetParLimits(kmu25,2.e-3,3.e-3);
  f.SetParLimits(ks25,+1.e-3,3.e-3);
  
  f.FixParameter(kmu25,1.5e-3);
  f.SetParameter(ks25,1.2e-3);
  
  
  f.SetParLimits(kmu50,4.e-3,6.e-3);
  f.SetParLimits(ks50,+1.e-3,4.e-3);
  
  f.FixParameter(kmu50,5.0e-3);
  f.SetParameter(ks50,2.5e-3);
  
  
  f.SetParLimits(kmuO,-4.e-4,4.e-4);
  f.SetParLimits(ksO,+3.e-3,1.e-2);
  
  f.FixParameter(kmuO,0.e-3);
  f.SetParameter(ksO,4.5e-3);

  f.SetParLimits(knorm,1.e-2,1.e+3);
  f.SetParameter(knorm,.6);

  std::cout << "All parameters fixed...\n";

  TCanvas *c1 = new TCanvas("c1","",800,600);
  c1->cd();

  gStyle->SetOptFit(1111);

  uRes->GetXaxis()->SetRangeUser(xmin,xmax);
  uRes->Draw();
  uRes->Fit("f","Rv");


  TString uname(Form("plots/BAD_ResU_Run%d.png",runNumber));

  c1->SaveAs(uname.Data());

  vRes->GetXaxis()->SetRangeUser(xmin,xmax);
  vRes->Draw();
  vRes->Fit("f","Rv");


  TString vname(Form("plots/BAD_ResV_Run%d.png",runNumber));

  c1->SaveAs(vname.Data());


  

}