#include <math.h>

TF1 *fitFunc;
Int_t System=0; 
Int_t PART;
TF1 *fplf;  
//             0  1  2  3   4  5  6   7   8   9  10  11
//             P, D, T,3He,4He,Li,Be, B,  C,  N  O    F
Int_t Mass[] ={1.,2.,3.,3.,4.,7.,10.,11.,13.,14.,16.,18.};
Int_t Z[]    ={1,1,1,2,2,3,4,5,6,7,8,9};
Int_t Sys[]  ={86,86,78,78};
const Char_t *reaction[] = {"86Kr64Ni","86Kr58Ni","78Kr58Ni","78Kr64Ni"};
const Char_t *title[] = {"^{86}Kr+^{64}Ni","^{86}Kr+^{58}Ni","^{78}Kr+^{58}Ni","^{78}Kr+^{64}Ni"};

Int_t verbose;

void FitT(Int_t partId) {
  PART = partId;

  TString reactadr = Form("c.root");
  //TString reactadr = Form("/data/sjygroup/sjy11/wuenschel/040607/PhysTapes/SourceDef/%s_Z3034_Vfcut_Qcut_N.root",reaction[System]);

  TFile f(reactadr.Data());
  if(!f.IsOpen()){
    printf("file not opened: %s\n",reactadr.Data());
    return;
  }
  const Char_t *part[] = {"p","d","t","^{3}He","^{4}He","Li","Be","B","C","N","O","F"};
  const Char_t *matrpart[] = {"P0","P1","P2","P3","P4","Z3","Z4","Z5","Z6","Z7","Z8","Z9"};
  TCanvas *c = new TCanvas("c"," c",1100,700);
  c->Divide(4,3);
  TGraphErrors *Mpar = new TGraphErrors(20);
  TGraphErrors *Bpar = new TGraphErrors(20);
  TGraphErrors *Tpar = new TGraphErrors(20);     
  TGraphErrors *Ppar = new TGraphErrors(20);  
  Double_t ref =0.;     
       
  fplf = new TF1("fplf",funx_plf_source,0.,120.,4.);
  fplf->SetParNames("M","T","B","P");

  TString hfname;
  //hfname=Form("/home/sjygroup/wuenschel/040607/dMovingSource/StartParameters_%s_Z3.h",reaction[System]);
  hfname=Form("StartParameters_%s_%s.h",reaction[System],matrpart[partId]);
  gROOT->LoadMacro(hfname.Data());
  
  for(Int_t i=0; i<4; i++){
     printf("%s %lf %lf %lf\n",parName[i]->Data(),vstart[i],parMin[i],parMax[i]);
     }
 
  f.Cd("Edist");

  Double_t last=0.; 

  // for(Int_t ii=2; ii<13; ii++){
  //for(Int_t ii=1; ii<8; ii++){
  for(Int_t ii=1; ii<2; ii++){
     //   c->cd(ii-1);
     c->cd(ii);
     TH1F    *hd;

     TString hName = Form("eDist_E*%d_%s",ii,matrpart[partId]); 
     printf("%s\n",hName.Data());
       
     hd = (TH1F*)gROOT->FindObject(hName.Data());
     hd->SetDirectory(0);

     hd->Sumw2();
//      for(Int_t bin =0; bin<hd->GetNbinsX(); bin++){
//         Double_t yield = hd->GetBinContent(bin);
//         if(0.1*yield > hd->GetBinError(bin)){
// 	 hd->SetBinError(bin,0.1*yield);
//            }
//         }
         
       if(hd->Integral()<100.) continue;
       hd->Scale(1./hd->Integral());
       
       hd->SetMarkerStyle(20);
       hd->SetMarkerColor(kBlack);
       hd->SetMarkerSize(0.5);
       hd->SetLineColor(kBlack);
       hd->Draw("p");
       //gPad->SetLogy(kTRUE);
       
       
       //Now, do the fit.
       //First, set up fit function
       for(Int_t j=0;j<4;j++) {
          fplf->SetParameter(j,vstart[j]);  //Initial parameters
          }
       
       hd->Fit(fplf,"","",1,100.);
//        Double_t matrix[4][4];
//        gMinuit->mnemat(&matrix[0][0],4);

//        for(Int_t j=0; j<4; j++){
//          printf("%lf %lf %lf %lf\n",matrix[j][0],matrix[j][1],matrix[j][2],matrix[j][3]);
//        }

       //Set parameters in functions
       Double_t par[5],parError[5];
       for(i=0;i<4;i++) {
          par[i] = fplf->GetParameter(i);
	parError[i] = fplf->GetParError(i);
          }

       // ref = 2.477;
       if(ii==1) ref = par[2];
       Double_t En = Double_t(ii)+0.5;
       //Double_t En = Double_t(ii)+0.25;
       //   Double_t factor = (Double_t(Sys[System])-2.*Mass[partId])/(Double_t(Sys[System])-Mass[partId]);
       Double_t factor = 1.;
       Mpar->SetPoint(ii,En,par[0]);
       Mpar->SetPointError(ii,0.,parError[0]);
       Tpar->SetPoint(ii,En,par[1]/factor);
       Tpar->SetPointError(ii,0.,parError[1]);
       Double_t density;
       if(ref>0.) density = par[2]/ref;
       //Bpar->SetPoint(ii,Double_t(ii),par[2]);
       Bpar->SetPoint(ii,En,density);
       Bpar->SetPointError(ii,0.,parError[2]);
       Ppar->SetPoint(ii,En,par[3]);
       Ppar->SetPointError(ii,0.,parError[3]);

       TF1 *fplf1 = fplf->Clone();
       fplf1->SetLineColor(kBlue);
       
       fplf1->Draw("same");
       fplf1->GetXaxis()->SetRangeUser(0,120);

       printf("\n\n end of %d ref %lf\n\n",ii,ref);
     
     }

     TCanvas *c2 = new TCanvas("par","par",600,600);
     c2->cd();
     Tpar->SetMarkerStyle(23);
     Tpar->SetMarkerColor(1);
     Tpar->Draw("AP");
     Mpar->SetMarkerColor(2);
     Mpar->SetMarkerStyle(23);
     Mpar->Draw("sameP");
     Bpar->SetMarkerColor(3);
     Bpar->SetMarkerStyle(23);
     Bpar->Draw("sameP"); 
     Ppar->SetMarkerColor(4);
     Ppar->SetMarkerStyle(23);
     Ppar->Draw("sameP");
     TLegend *leg = new TLegend(.14,.7,.2,.85);
     leg->AddEntry(Mpar,"M","P");
     leg->AddEntry(Bpar,"B","P");
     leg->AddEntry(Tpar,"T","P");
     leg->AddEntry(Ppar,"P","P");
     leg->SetFillColor(0);
     leg->SetBorderSize(1);
     leg->SetTextFont(102);
     leg->SetTextSize(0.02);
     leg->Draw("same");

}

Double_t funx_plf_source(Double_t *x,Double_t *par) {
  Float_t PI = TMath::Pi();
   
  ///////////////////////////////////////////////////////////////////////////
    // PLF source
    Float_t Mc  = par[0];
    Float_t Tc  = par[1];
    Float_t Bcc = par[2];
    //Float_t Pc  = 1.5;
    Float_t Pc  = par[3];

    ///////////////////////////////////////////////////////////////////////////
     
      Double_t diff = x[0] - Bcc;
      Double_t resc = Mc*((2*diff - Pc)*TMath::Exp(-diff/Tc)*TMath::Erfc((Pc-2*diff)/(2*TMath::Sqrt(Pc*Tc)))+2*TMath::Sqrt((Pc*Tc)/PI)*TMath::Exp(-(Pc*Pc+4*(diff*diff))/(4*Pc*Tc)));
      //Double_t resc = Mc*diff*TMath::Exp(-diff/Tc);

        // Total
        Float_t res = resc;
        
        return res;

}