TCanvas *c;
double par[6];  
double *epar;  
float mpvInc; // this is ugly


pl3() {
  gStyle->SetStatW(0.22);
  gStyle->SetStatH(0.22);
  f=new TFile("june18.hist.root");

  TH1F *h1=new TH1F("mpv","MPV from Pres-1, MIP gated; MPV of ADC-ped",25,-5,50); 
  TH1F *h2=new TH1F("mpvE","relative error of MPV from Pres-1, MIP gated; err(MPV)/MPV ",25,0,0.6); 

  TH1F *h3=new TH1F("mpvS","relative width of L-peak, from Pres-1, MIP gated; sigma/MPV",25,0.,1.); 


 TH2F *h4=new TH2F("mpv2","MPV from Pres-1; gated w/ MIP ; inclusive spectrum;",25,0,50,25,0,25); 
 


  char core[100];
  int eta;
  char cT='P';
  char sub='C';

  for(sub='A';sub<='E';sub++)
  for(eta=1;eta<=12;eta++) {
    sprintf(core,"05%c%c%02d",cT,sub,eta);
    TString coreT=core;
    ha=(TH1F*)f->Get("a"+coreT);
    h=(TH1F*)f->Get("j"+coreT);
    c=new TCanvas("aa","aa",400,400);;
    plotOne(ha,h);
    c->Print(coreT+".ps");
    c->Print(coreT+".gif");
    float mpv=par[4];
    int mpvEr=epar[4];

    h1->Fill(mpv);
    if(mpv>0) {
      h2->Fill(mpvEr/mpv);
      h3->Fill(par[5]/mpv);
    }
    h4->Fill(mpv,mpvInc);

  }
  gStyle->SetOptStat(1111111);

  c=new TCanvas("res","res",600,600);
  c->Divide(2,2);
  c->cd(1); h1->Draw();
  c->cd(2); h2->Draw();
  c->cd(3); h3->Draw();
  c->cd(4); h4->Draw("box");

}

//=================================
plotOne(TH1F *ha, TH1F *h) {


  float xMax=120;

  h->SetAxisRange(5,40);
  float sum=h->Integral();
  printf("sum=%f\n",sum);
  if(sum>110) 
    h->Rebin(2);
  else if(sum>70)
    h->Rebin(3);
  else 
    h->Rebin(4);

  ha->SetAxisRange(-10,xMax); 
  h->SetAxisRange(-10,xMax);


  c->Divide(1,2);

  c->cd(1);
  ha->Draw();
  // gPad->SetGrid();
  float yMax=1.5*fitGausLand(ha);
  mpvInc=par[4];
  if(yMax<300) yMax=1200;
  ha->SetMaximum(yMax);
 

  c->cd(2);
  h->SetMaximum(1.2*h->GetMaximum());
  fitGausLand(h);
  // gPad->SetGrid();

}


//---------------------------------------
float fitGausLand(TH1F *h) {

  float xMax=80.;

  gF=new TF1("gF","gaus",-2,2);
  gF->SetLineWidth(2);
  gF->SetLineColor(kBlue);
  gF->SetLineStyle(2);

  lF=new TF1("lF","landau",3,xMax);
  lF->SetLineWidth(2);
  lF->SetLineStyle(2);
  lF->SetLineColor(kRed);

  glF=new TF1("glF","gaus(0)+landau(3)",-8,xMax);
  glF->SetLineWidth(1);
  glF->SetLineColor(kGreen);
  glF->SetParNames("ampl-G","mean-G","sig-G","const-L","MPV-L","sig-L");


  h->Fit("gF","RQ");
  // return 1;
  h->Fit("lF","R0Q+");
  
  // copy starting point
  gF->GetParameters(par+0); par[0]/=2;
  lF->GetParameters(par+3);
  glF->SetParameters(par);

  h->Fit("glF","R");

  // retrieve components
  glF->GetParameters(par);
  epar= glF->GetParErrors();

  gF->SetParameters(par+0);
  lF->SetParameters(par+3);

  gF->SetRange(-5,5);
  gF->Draw("same");

  lF->SetRange(0,xMax);
  lF->Draw("same");

  float yMax=lF->GetMaximum();
  printf("max=%f\n",yMax);
  return yMax;
}