FitMass.C

void Fit(TH1F* hmass, string Energy);
TFile* ofile;
int rebin=1;

//////////////////////////////////////////////////////
void FitMass(){

  TFile* ifile = new TFile("histo_mass_all.root","read");
  //TFile* ifile = new TFile("histo_mass_part2.root","read");

  ofile = new TFile("histo_mass_fitted.root","recreate");

  for(int i=0; i<11; i++){
  //for(int i=2; i<3; i++){
    int Ex = i+6;
    if(i>4) rebin=2;
    else rebin=1;

    string Energy = to_string(Ex) + "MeV";
    TString histo_name = Form("hmass_%iMeV",i+6);
    TH1F* h1 = (TH1F*)ifile->FindObjectAny(histo_name);
    Fit(h1,Energy);
  }

  ofile->Close();
}

//////////////////////////////////////////////////////
void Fit(TH1F* hmass,string Energy){
  hmass->Rebin(rebin);
  hmass->Draw("E1");

  int Amin = 82;
  int Amax = 153;
  //int Amin = 95;
  //int Amax = 106;


  int A = Amin;

  int NumberOfA = Amax - Amin;
  Double_t para[3*NumberOfA];

  TF1* g[NumberOfA];
  double Integral[NumberOfA];
  double Integral_err[NumberOfA];
  double total_integral = 0;
  double Yield[NumberOfA];

  for(int i=0; i<NumberOfA; i++){
    g[i] = new TF1(Form("g%i",i),"gaus",A-0.4,A+0.4);   
    if(i==0){
      g[i]->SetParameter(0,100);
      g[i]->SetParameter(1,Amin);
      g[i]->SetParameter(2,0.3);
    }
    else if(i>0){
      g[i]->SetParameter(0,g[i-1]->GetParameter(0));
      g[i]->SetParameter(1,A);
      g[i]->SetParameter(2,g[i-1]->GetParameter(2));
    }
 
    g[i]->SetParLimits(1,A-0.1,A+0.1);
    g[i]->SetParLimits(2,0.25,0.37);

   
    hmass->Fit(g[i],"qBR");

    g[i]->Draw("lsame");
    g[i]->GetParameters(&para[3*i]);

    //cout << A << " " << g[i]->GetParameter(2) << endl;
    A++;
  }


  TString total_func = "gaus(0)";
  for(int i=1; i<NumberOfA; i++){
    TString gaus_func = Form("gaus(%i)",3*i);
    total_func += "+" + gaus_func;
  }

  TF1* total = new TF1("total",total_func,Amin-0.5, Amin+NumberOfA-0.5);
  total->SetParameters(para);
  A = Amin;
  for(int i=0; i<NumberOfA; i++){
    total->SetParameter(3*i,g[i]->GetParameter(0));
    total->SetParameter(3*i+1,g[i]->GetParameter(1));
    total->SetParameter(3*i+2,g[i]->GetParameter(2));
    total->SetParLimits(3*i+1,A-0.2,A+0.2);
    total->SetParLimits(3*i+2,g[i]->GetParameter(2)*0.85,g[i]->GetParameter(2)*1.01);
    A++;
  }
  total->SetLineColor(4);
  total->SetNpx(1000);

  hmass->Fit(total,"RBq");

  A = Amin;
  //TF1* one_gauss = new TF1("one_gaus","gaus");
  for(int i=0; i<NumberOfA; i++){
      
    double Amplitude = total->GetParameter(3*i);
    double mean = total->GetParameter(3*i+1);
    double sigma = total->GetParameter(3*i+2);
    double Amplitude_err = total->GetParError(3*i);
    double mean_err = total->GetParError(3*i+1);
    double sigma_err = total->GetParError(3*i+2);
 
    /*double Amplitude = g[i]->GetParameter(0);
    double mean = g[i]->GetParameter(1);
    double sigma = g[i]->GetParameter(2);
    double Amplitude_err = g[i]->GetParError(0);
    double mean_err = g[i]->GetParError(1);
    double sigma_err = g[i]->GetParError(2);*/
    
    //one_gauss->SetParameter(0,Amplitude);
    //one_gauss->SetParameter(1,mean);
    //one_gauss->SetParameter(2,sigma);
    //one_gauss->SetParError(0,Amplitude_err);
    //one_gauss->SetParError(1,mean_err);
    //one_gauss->SetParError(2,sigma_err);

    cout << "************************" << endl;
    cout << "A=  " << A << endl;
    //cout << "Amplitude= " << Amplitude << "+/-" << Amplitude_err << endl;
    //cout << "Mean= " << mean  << endl;
    cout << "Sigma= " << sigma << "+/-" << sigma_err << endl;

    

    Integral[i] = (12.5/rebin)*Amplitude*sigma*sqrt(2*TMath::Pi());
    Integral_err[i] = (12.5/rebin)*sqrt(2*TMath::Pi()*(pow(sigma*Amplitude_err,2) + pow(Amplitude*sigma_err,2)));
    //Integral_err[i] = sqrt(Integral[i]);
    

    total_integral += Integral[i];
    A++;
  }

  double histo_integral = hmass->Integral();
  hmass->Write();
  total->Write();

  string output_filename = "Mass_Yield_" + Energy + ".dat";
  ofstream ofile;
  ofile.open(output_filename.c_str());

  // Yield calculation //
  A = Amin;
  for(int i=0; i<NumberOfA; i++){
    //Yield[i] = Integral[i]/total_integral*200;
    Yield[i] = Integral[i]/histo_integral*200;
    //double yield_err = Integral_err[i]/total_integral*200;
    double yield_err = Integral_err[i]/histo_integral*200;
    ofile << A << " " << Yield[i] << " " << yield_err << endl;
    if(A==99) cout << "A=99 -> Y= " << Yield[i] << endl;
    if(A==147) cout << "A=147 -> Y= " << Yield[i] << endl;
    A++;
  }
  ofile.close();


}