const unsigned int n_bins(100); //const double m_min(0.1399), m_max(0.16); const double m_min(0.1399), m_max(0.159); const double bin_width = (m_max-m_min)/((double)n_bins); double gaus_func(double m, double mu, double sigma) { if (mm_max) return 0.; double norm = TMath::Sqrt(TMath::Pi()/2.)*sigma*(TMath::Erf((mu-m_min)/(TMath::Sqrt2()*sigma))-TMath::Erf((mu-m_max)/(TMath::Sqrt2()*sigma))); if (norm==0.) return 0.; return TMath::Gaus(m,mu,sigma,0)/norm; } double sgn_func(double *x, double *p) { double m = x[0]; if (mm_max) return 0.; double Nsig = p[0]; double mu = p[1]; double sigma= p[2]; return Nsig*(gaus_func(m,mu,sigma))*bin_width; } double bkg_func(double *x, double *p) { double m = x[0]; if (mm_max) return 0.; double Nbkg = p[0]; double alpha = p[1]; double beta = p[2]; double m0 = p[3]; double norm = ((2./3.) * (pow((m_max-m0),3./2.)) + (2./5.) * alpha* (pow((m_max-m0),5./2.)) + (2./7.) * beta *(pow((m_max-m0),7./2.))) - ((2./3.) * (pow((m_min-m0),3./2.)) + (2./5.) * alpha* (pow((m_min-m0),5./2.)) + (2./7.) * beta *(pow((m_min-m0),7./2.))); if (norm==0.) return 0.; return Nbkg*((pow((m-m0),1./2.) + alpha*pow((m-m0),3./2.) + beta*pow((m-m0),5./2.))/norm)*bin_width; } double tot_func(double *x, double *p) { if (x[0]m_max) return 0.; return sgn_func(x,p)+ bkg_func(x,&p[3]); } void fit_mc() { //TFile * f = new TFile("check_here_for_the_plots_in_charmWG.root"); TFile * f = new TFile("histo_ksks_proc11_bucket_9to12_unscaled.root"); //TFile * f = new TFile("histo_ksks_proc11_unscaled.root"); //TFile * f = new TFile("histo_ksks_proc11_scaled.root"); //TFile * f = new TFile("histo_ksks_bucket_9to12_unscaled.root"); //TFile * f = new TFile("histo_ksks_bucket_9to12_scaled.root"); //TFile * f = new TFile("histo_ksks_proc11_bucket_9to12_scaled.root"); TH1F *h1 = new TH1F("h1","histogram",100, 0.139, 0.159); // h1 = (TH1F*)f-> Get("DeltaM_fde"); h1 = (TH1F*)f-> Get("DeltaM4"); gStyle->SetOptFit(1111); gStyle->SetHistMinimumZero(kTRUE); // Fit TF1 *tot = new TF1("tot",tot_func,m_min,m_max,7); tot->SetNpx(n_bins*3); //tot->SetParameters(2e3,0.145,0.0006,1e2,9,20,0.1399);//ft >0.001 to 0.003 tot->SetParameters(2e3,0.145,0.0006,1e2,9,20,0.1399);//fte > 4 tot->SetParNames("N_{sig}","#mu","#sigma","N_{bkg}","#alpha","#beta","m0"); tot->SetParLimits(1,m_min,m_max);//mu1 tot->SetParLimits(2,0.,10.);//sigma1 tot->FixParameter(6,0.1399); //m0 tot->FixParameter(5,0.); //beta //h1->Fit(tot,"R"); h1->Fit(tot,"L"); h1->SetMarkerStyle(20); h1->Draw("E"); TF1 *sgn = new TF1("sgn",sgn_func,m_min,m_max,3); sgn->SetLineColor(kGreen); sgn->SetLineStyle(2); sgn->SetParameters(tot->GetParameter(0),tot->GetParameter(1),tot->GetParameter(2)); //sgn->SetParameters(tot->GetParameters()); sgn->Draw("same"); TF1 *bkg = new TF1("bkg",bkg_func,0.139,0.159,4); bkg->SetLineColor(kBlue); bkg->SetLineStyle(2); bkg->SetParameters(tot->GetParameter(3),tot->GetParameter(4),tot->GetParameter(5),tot->GetParameter(6)); bkg->Draw("same"); double w=h1->GetBinWidth(1); double integral0 = sgn->Integral(0.1445,0.1465)/w; double integral1 = bkg->Integral(0.1445,0.1465)/w; cout<<"Integral of sgn in [0.1445,0.1465]= "<< integral0 <