#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TF1.h"
#include "TMinuit.h"
#include "TMath.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TLatex.h"


static const double lam0 = 0.4551599;
static const double LAMS = 0.000688;
static const double LAMT = 0.000012;
static const double lam_of = 2.5;
static const double lam_pf = 0.0074;
static const double lam_op = 0.041;
static const double LAM_OM = 1.009 *(0.75*LAMS + 0.25*LAMT) ;
static const double LAM_PM = 1.143 *(0.25*LAMS + 0.75*LAMT) ;


/* *********************************************************************
********************************************************************* */
double n1_population_function(double *var, double *par)
{
  double time = var[0];
  double phi = par[0];
  double lam1 = LAMS + phi*lam_of + phi*lam_pf;
  double n1 = exp(-lam1*time);
  return n1;
}


/* *********************************************************************
********************************************************************* */
double n2_population_function(double *var, double *par)
{
  double time = var[0];
  double phi = par[0];
  double lam1 = LAMS + phi*lam_of + phi*lam_pf;
  double lam2 = LAM_OM + lam_op;
  double n2 = phi*lam_of/(lam1-lam2) * 
    (exp(-lam2*time) - exp(-lam1*time));
  return n2;
}


/* *********************************************************************
********************************************************************* */
double n3_population_function(double *var, double *par)
{
  double time = var[0];
  double phi = par[0];
  double lam1 = LAMS + phi*lam_of + phi*lam_pf;
  double lam2 = LAM_OM + lam_op;
  double lam3 = LAM_PM;
  double n3a = phi*lam_of*lam_op / (lam1-lam2) * 
    (((exp(-lam3*time) - exp(-lam2*time))/(lam2-lam3)) 
     - ((exp(-lam3*time) - exp(-lam1*time))/(lam1-lam3)));
  double n3b = phi*lam_pf / (lam1-lam3) * (exp(-lam3*time) - exp(-lam1*time));
  double n3 = n3a + n3b;
  return n3;
}


/* *********************************************************************
********************************************************************* */
void plot_populations()
{
  gROOT->SetStyle("Plain");
  TCanvas *c1 = new TCanvas("canvas", "canvas", 800, 400);
  c1->Divide(3,1);
  TF1* n[4][4];

  for (int index=1; index<=3; index++) {

    double phi;
    switch(index) {
    case 1: phi = 1.0; break;
    case 2: phi = 0.1; break;
    case 3: phi = 0.01; break;
    }

    char n1name[20], n2name[20], n3name[20];
    sprintf(n1name, "n1_%d", index);
    sprintf(n2name, "n2_%d", index);
    sprintf(n3name, "n3_%d", index);
    n[1][index] = new TF1(n1name, n1_population_function, 0, 25, 1);
    n[2][index] = new TF1(n2name, n2_population_function, 0, 25, 1);
    n[3][index] = new TF1(n3name, n3_population_function, 0, 25, 1);
    for (int nlabel=1; nlabel<=3; nlabel++) {
      n[nlabel][index]->SetNpx(1000);
      n[nlabel][index]->FixParameter(0,phi);
    }

    c1->cd(index);
    n[1][index]->UseCurrentStyle();
    n[1][index]->SetMinimum(0);
    n[1][index]->SetMaximum(1);
    n[1][index]->GetXaxis()->SetLabelOffset(-0.07);
    n[1][index]->GetXaxis()->SetTicks("-");
    n[1][index]->GetXaxis()->SetNdivisions(006);
    n[1][index]->GetYaxis()->SetLabelOffset(-0.04);
    n[1][index]->GetYaxis()->SetTicks("+");
    n[1][index]->GetYaxis()->SetNdivisions(505);


    gStyle->SetLineStyleString(12,"40 30");
    gStyle->SetLineStyleString(13,"5 15");

    n[1][index]->SetLineStyle(1);
    n[2][index]->SetLineStyle(12);
    n[3][index]->SetLineStyle(13);

    n[1][index]->SetLineWidth(1);
    n[2][index]->SetLineWidth(1);
    n[3][index]->SetLineWidth(1);

    n[1][index]->Draw("C");
    n[2][index]->Draw("CSAME");
    n[3][index]->Draw("CSAME");

  }
}


/* ****************************************************************** */