#include <fstream>
void fit_peaks(){
  TH1F *h2GammaMassECAL=(TH1F*)_file0->FindObjectAny("h2GammaMassECAL");
  h2GammaMassECAL->GetXaxis()->SetRangeUser(0.08,0.2);
  TF1 *f2=new TF1("f2","gaus(0)+gaus(3)+expo(6)",0,0.4);
  double max=h2GammaMassECAL->GetMaximum();
  if (max>100){
    f2->SetParameters(100.,0.13,0.0035,10.,0.13,0.007,6.3,-25.5);
  }
  else {
    f2->SetParameters(50.,0.13,0.004);	
  }
  h2GammaMassECAL->Fit("f2","q");
  if (f2->GetParameter(0)>10.){
    double sigma1=1000.*f2->GetParameter(2);
    double sigma2=1000.*f2->GetParameter(5);
    double N1=f2->GetParameter(0);
    double N2=f2->GetParameter(3);
    double avg_sigma=(sigma1*N1+sigma2*N2)/(N1+N2);
  
    ofstream out("m.dat",ios::app);
    out << 1000.*f2->GetParameter(1) << " " <<  avg_sigma;
  
    TH1F *h2GammaMassFCAL=(TH1F*)_file0->FindObjectAny("h2GammaMassFCAL");
    h2GammaMassFCAL->GetXaxis()->SetRangeUser(0.08,0.2);
    f2->SetParameters(500.,0.13,0.007,100.,0.13,0.02,4.77,-7.6);
    if (max>100){
      f2->SetParameters(500.,0.13,0.007,100.,0.13,0.02,4.77,-7.6);
    }
    else {
      f2->SetParameters(50.,0.13,0.007,0.,0.,0.,0.,0.);	
    }
    h2GammaMassFCAL->Fit("f2","q");
    sigma1=1000.*f2->GetParameter(2);
    sigma2=1000.*f2->GetParameter(5);
    N1=f2->GetParameter(0);
    N2=f2->GetParameter(3);
    avg_sigma=(sigma1*N1+sigma2*N2)/(N1+N2);
    
    out << " " <<1000.*f2->GetParameter(1) << " " << avg_sigma<< endl;
  }
}