class MyFitFun1 { public: double Evaluate(double *xPtr, double *parPtr) { double fitVal = 0.0; double xVal = *xPtr; double gMag = parPtr[0]; double center = parPtr[1]; double sigma = parPtr[2]; double gFun = TMath::Gaus(xVal,center,sigma); //Background fit parameters double p0Val = parPtr[3]; double p1Val = parPtr[4]; double p2Val = parPtr[5]; //Background function double bkg = p0Val + p1Val*xVal + p2Val*xVal*xVal; //Fit value fitVal = gMag*gFun + bkg; return fitVal; } }; class MyFitFun2 { public: double Evaluate(double *xPtr, double *parPtr) { double fitVal = 0.0; double xVal = *xPtr; double gMag = parPtr[0]; double center = parPtr[1]; double sigma = parPtr[2]; double gFun = TMath::Gaus(xVal,center,sigma); //Background fit parameters double p0Val = parPtr[3]; double p1Val = parPtr[4]; double p2Val = parPtr[5]; //Background function double bkg = fabs(p0Val + p1Val*xVal + p2Val*xVal*xVal); //a0(980) double mag2 = parPtr[6]; double center2 = parPtr[7]; double width2 = parPtr[8]; //double gFun2 = mag2*TMath::Gaus(xVal,center2,width); //double bwFun = mag2/(pow(xVal*xVal - center2*center2,2) + pow(center2*width2,2)); double voigtFun = TMath::Voigt(xVal-center2,0.006,width2); //double sFun = 1.0/(exp((-xVal + 1.0)*200) + 1); double sCenter = parPtr[9]; double sWidth = parPtr[10]; double sFun = 1.0/(exp((-xVal + sCenter)*sWidth) + 1); //fitVal = sFun*gMag*gFun + sFun*bkg + sFun*bwFun; fitVal = sFun*gMag*gFun + sFun*bkg + sFun*voigtFun; return fitVal; } }; void fitPhiPeaks(int iShow){ gROOT->Reset(); gROOT->SetStyle("Plain"); gStyle->SetOptStat(0); char hId[120]; TFile *hFile = new TFile("./all_kpkmpi0.root");//CHANGE ME TH2D *hMassXY; hFile->GetObject("hKpKmVsKpKmPi0CutD4",hMassXY);//CHANGE ME //TH2D *hMassXY; hFile->GetObject("hKpKmVsKpKmPi0",hMassXY);//CHANGE ME MyFitFun1 * fptr1 = new MyFitFun1(); int nPar = 6; TF1 *f1 = new TF1("f1",fptr1,&MyFitFun1::Evaluate,0.0,4.0,nPar, "MyFitFun1","Evaluate"); nPar = 11; MyFitFun2 * fptr2 = new MyFitFun2(); TF1 *f2 = new TF1("f2",fptr2,&MyFitFun2::Evaluate,0.0,4.0,nPar, "MyFitFun2","Evaluate"); TF1 *f0 = new TF1("f0","gaus"); f1->SetLineWidth(2); f1->SetLineColor(4); int nBinsX = hMassXY->GetNbinsX(); double xMin = hMassXY->GetXaxis()->GetXmin(); double xMax = hMassXY->GetXaxis()->GetXmax(); TH1D *hMassX = new TH1D("hMassX","",nBinsX,xMin,xMax); TH1D *hMassX1 = new TH1D("hMassX1","",nBinsX,xMin,xMax); TH1D *hMassX2 = new TH1D("hMassX2","",nBinsX,xMin,xMax); TH1D *hMassXB = new TH1D("hMassXB","",nBinsX,xMin,xMax); TH1D *hMassXP0 = new TH1D("hMassXP0","",nBinsX,xMin,xMax); TH1D *hMassXP1 = new TH1D("hMassXP1","",nBinsX,xMin,xMax); TH1D *hMassXP2 = new TH1D("hMassXP2","",nBinsX,xMin,xMax); TH1D *hMassXP3 = new TH1D("hMassXP3","",nBinsX,xMin,xMax); TH1D *hMassXP4 = new TH1D("hMassXP4","",nBinsX,xMin,xMax); TH1D *hMassXP5 = new TH1D("hMassXP5","",nBinsX,xMin,xMax); TH1D *hMassXP6 = new TH1D("hMassXP6","",nBinsX,xMin,xMax); TH1D *hMassXP7 = new TH1D("hMassXP7","",nBinsX,xMin,xMax); TH1D *hMassXP8 = new TH1D("hMassXP8","",nBinsX,xMin,xMax); TH1D *hMassXP9 = new TH1D("hMassXP9","",nBinsX,xMin,xMax); TH1D *hMassXP10 = new TH1D("hMassXP10","",nBinsX,xMin,xMax); int nBinsY = hMassXY->GetNbinsY(); double yMin = hMassXY->GetYaxis()->GetXmin(); double yMax = hMassXY->GetYaxis()->GetXmax(); TH1D *hMassY[nBinsX]; TH1D *hMassYSignal[nBinsX]; TH1D *hMassYSignal1[nBinsX]; TH1D *hMassYSignal2[nBinsX]; TH1D *hMassYBack[nBinsX]; TH1D *hMassYBack1[nBinsX]; TH1D *hMassYBack2[nBinsX]; int lBin = 1; int uBin = nBinsX; if (iShow > 0){ lBin = iShow; uBin = iShow; } for (int iBinX = lBin; iBinX <= uBin; iBinX++) { sprintf(hId,"hMassY_%d",iBinX); hMassY[iBinX-1] = (TH1D*) hMassXY->ProjectionY(hId,iBinX,iBinX,""); sprintf(hId,"hMassYSignal_%d",iBinX); hMassYSignal[iBinX-1] = new TH1D(hId,"",nBinsY,yMin,yMax); sprintf(hId,"hMassYSignal1_%d",iBinX); hMassYSignal1[iBinX-1] = new TH1D(hId,"",nBinsY,yMin,yMax); sprintf(hId,"hMassYSignal2_%d",iBinX); hMassYSignal2[iBinX-1] = new TH1D(hId,"",nBinsY,yMin,yMax); sprintf(hId,"hMassYBack_%d",iBinX); hMassYBack[iBinX-1] = new TH1D(hId,"",nBinsY,yMin,yMax); sprintf(hId,"hMassYBack1_%d",iBinX); hMassYBack1[iBinX-1] = new TH1D(hId,"",nBinsY,yMin,yMax); sprintf(hId,"hMassYBack2_%d",iBinX); hMassYBack2[iBinX-1] = new TH1D(hId,"",nBinsY,yMin,yMax); hMassY[iBinX-1]->GetXaxis()->SetRangeUser(0.96,1.2); double maxVal = hMassY[iBinX-1]->GetMaximum(); f0->SetParameter(0,maxVal); f0->SetParameter(1,1.02); f0->SetParameter(2,0.005); hMassY[iBinX-1]->Fit("f0","B","R",1.0,1.03); double par0 = f0->GetParameter(0); double par1 = f0->GetParameter(1); double par2 = f0->GetParameter(2); f1->SetParameter(0,par0); f1->FixParameter(2,7.12905e-03); f1->FixParameter(1,par1); f1->SetParameter(3,0); f1->SetParameter(4,0); f1->SetParameter(5,0); f1->SetParLimits(2,0.005,0.01); double fitHigh = 1.08; if (iBinX <= 66) fitHigh = 1.02; if (iBinX == 67) fitHigh = 1.03; if (iBinX == 68) fitHigh = 1.04; if (iBinX == 68 || iBinX == 69){ f1->FixParameter(4,0); f1->FixParameter(5,0); } f1->FixParameter(4,0); f1->FixParameter(5,0); hMassY[iBinX-1]->Fit("f1","B","R",0.98,fitHigh); double f1Par0 = f1->GetParameter(0); double f1Par1 = f1->GetParameter(1); double f1Par2 = f1->GetParameter(2); double f1Par3 = f1->GetParameter(3); double f1Par4 = f1->GetParameter(4); double f1Par5 = f1->GetParameter(5); cout<<"f1Par1 = "<SetParameter(0,f1Par0); //f2->SetParameter(1,f1Par1); f2->SetParameter(1,1.02); f2->FixParameter(2,7.12905e-03); f2->SetParameter(3,f1Par3); f2->SetParameter(4,f1Par4); //f2->SetParameter(5,f1Par5); //f2->FixParameter(3,0); //f2->FixParameter(4,0); f2->FixParameter(5,0); f2->SetParameter(6,0); //f2->SetParLimits(6,0.0,1000000000); f2->FixParameter(7,0.980); //f2->FixParameter(7,1.0); f2->FixParameter(8,0.075); //f2->SetParLimits(8,0.01,0.80); f2->SetParameter(9,0.99); f2->SetParameter(10,150); //SetParLimits f2->SetParLimits(0,0.0,400); f2->SetParLimits(1,1.0,1.03); //f2->SetParLimits(3,0.0,400); f2->SetParLimits(6,0.0,8); //f2->SetParLimits(7,0.97,1.03); //f2->SetParLimits(8,0.0,0.2); f2->SetParLimits(9,0.97,1.02); //f2->SetParLimits(9,0.986,0.987); f2->SetParLimits(10,100,500); //if (iBinX >= 96) f2->FixParameter(6,0); if (iBinX >= 0) f2->FixParameter(6,0); hMassY[iBinX-1]->Fit("f2","","R",0.98,fitHigh); double f2par0 = f2->GetParameter(0); double f2par1 = f2->GetParameter(1); double f2par2 = f2->GetParameter(2); double f2par3 = f2->GetParameter(3); double f2par4 = f2->GetParameter(4); double f2par5 = f2->GetParameter(5); double f2par6 = f2->GetParameter(6); double f2par7 = f2->GetParameter(7); double f2par8 = f2->GetParameter(8); double f2par9 = f2->GetParameter(9); double f2par10 = f2->GetParameter(10); //f2->FixParameter(0,f2par0); //f2->FixParameter(1,f2par1); //f2->FixParameter(2,f2par2); //f2->FixParameter(3,f2par3); //f2->FixParameter(4,f2par4); //f2->FixParameter(5,f2par5); //f2->ReleaseParameter(6); //f2->ReleaseParameter(7); //f2->ReleaseParameter(8); //f2->FixParameter(7,f2par7); //f2->FixParameter(8,f2par8); //f2->SetParLimits(8,0.05,0.2); //f2->FixParameter(9,f2par9); //f2->FixParameter(10,f2par10); //hMassY[iBinX-1]->Fit("f2","B","R",0.98,fitHigh); double amp1 = f2->GetParameter(0); double amp1Err = f2->GetParError(0); double amp2 = f2->GetParameter(3); double amp2Err = f2->GetParError(3); hMassXP0->SetBinContent(iBinX,amp1); hMassXP0->SetBinError(iBinX,amp1Err); hMassXP1->SetBinContent(iBinX,f2->GetParameter(1)); hMassXP1->SetBinError(iBinX,f2->GetParError(1)); hMassXP2->SetBinContent(iBinX,f2->GetParameter(2)); hMassXP2->SetBinError(iBinX,f2->GetParError(2)); hMassXP3->SetBinContent(iBinX,amp2); hMassXP3->SetBinError(iBinX,amp2Err); hMassXP4->SetBinContent(iBinX,f2->GetParameter(4)); hMassXP4->SetBinError(iBinX,f2->GetParError(4)); hMassXP5->SetBinContent(iBinX,f2->GetParameter(5)); hMassXP5->SetBinError(iBinX,f2->GetParError(5)); hMassXP6->SetBinContent(iBinX,f2->GetParameter(6)); hMassXP6->SetBinError(iBinX,f2->GetParError(6)); hMassXP7->SetBinContent(iBinX,f2->GetParameter(7)); hMassXP7->SetBinError(iBinX,f2->GetParError(7)); hMassXP8->SetBinContent(iBinX,f2->GetParameter(8)); hMassXP8->SetBinError(iBinX,f2->GetParError(8)); hMassXP9->SetBinContent(iBinX,f2->GetParameter(9)); hMassXP9->SetBinError(iBinX,f2->GetParError(9)); hMassXP10->SetBinContent(iBinX,f2->GetParameter(10)); hMassXP10->SetBinError(iBinX,f2->GetParError(10)); double center = f1->GetParameter(1); double sigma = f1->GetParameter(2); double lowVal = center - 3*sigma; double highVal = center + 3*sigma; lowVal = 0.98; highVal = fitHigh; int binLow = hMassY[iBinX-1]->GetXaxis()->FindBin(lowVal); int binHigh = hMassY[iBinX-1]->GetXaxis()->FindBin(highVal); f2par0 = f2->GetParameter(0); f2par1 = f2->GetParameter(1); f2par2 = f2->GetParameter(2); f2par3 = f2->GetParameter(3); f2par4 = f2->GetParameter(4); f2par5 = f2->GetParameter(5); f2par6 = f2->GetParameter(6); f2par7 = f2->GetParameter(7); f2par8 = f2->GetParameter(8); f2->SetParameter(0,0.0); f2->SetParameter(3,0.0); f2->SetParameter(4,0.0); f2->SetParameter(5,0.0); hMassYSignal1[iBinX-1]->Eval(f2);//a0(980) f2->SetParameter(0,f2par0); f2->SetParameter(6,0); hMassYSignal2[iBinX-1]->Eval(f2);//phi f2->SetParameter(0,0.0); f2->SetParameter(6,0.0); f2->SetParameter(3,f2par3); f2->SetParameter(4,f2par4); f2->SetParameter(5,f2par5); hMassYBack[iBinX-1]->Eval(f2); hMassYSignal[iBinX-1]->Add(hMassY[iBinX-1],hMassYBack[iBinX-1],1.0,-1.0); hMassYBack[iBinX-1]->SetLineWidth(2); hMassYBack[iBinX-1]->SetLineColor(2); double signalVal = hMassYSignal[iBinX-1]->Integral(binLow,binHigh); double allVal = hMassY[iBinX-1]->Integral(binLow,binHigh); double backVal = hMassYBack[iBinX-1]->Integral(binLow,binHigh); double varSignal = allVal + backVal; double sigmaSignal = sqrt(varSignal); double signalVal1 = hMassYSignal1[iBinX-1]->Integral(binLow,binHigh); double signalVal2 = hMassYSignal2[iBinX-1]->Integral(binLow,binHigh); if (varSignal > 0) sigmaSignal = sqrt(varSignal); hMassYSignal[iBinX-1]->GetXaxis()->SetRangeUser(lowVal,highVal); hMassYSignal[iBinX-1]->SetFillStyle(3001); hMassYSignal[iBinX-1]->SetFillColor(3); hMassYSignal[iBinX-1]->SetLineColor(3); hMassYSignal1[iBinX-1]->SetLineColor(4); hMassYSignal2[iBinX-1]->SetLineColor(7); hMassYSignal1[iBinX-1]->SetLineWidth(2); hMassYSignal2[iBinX-1]->SetLineWidth(2); hMassX->SetBinContent(iBinX,signalVal); hMassX->SetBinError(iBinX,sigmaSignal); hMassX1->SetBinContent(iBinX,signalVal1); double frac1Err = 0; if (amp1 !=0) frac1Err = amp1Err/amp1; double signalErr1 = signalVal1*frac1Err; hMassX1->SetBinError(iBinX,signalErr1); hMassX2->SetBinContent(iBinX,signalVal2); double frac2Err = 0; if (amp2 !=0) frac2Err = amp2Err/amp2; double signalErr2 = signalVal2*frac2Err; hMassX2->SetBinError(iBinX,signalErr2); hMassXB->SetBinContent(iBinX,backVal); hMassXB->SetBinError(iBinX,sqrt(backVal)); cout<<"signal = "<SetMinimum(-50); hMassY[iBinX-1]->GetXaxis()->SetRangeUser(lowVal,highVal); hMassY[iBinX-1]->Draw("e"); hMassYBack[iBinX-1]->Draw("csame"); hMassYSignal[iBinX-1]->Draw("same"); hMassYSignal1[iBinX-1]->Draw("csame"); hMassYSignal2[iBinX-1]->Draw("csame"); //hMassX1->Draw("e1"); double binCenter = hMassX->GetBinCenter(iBinX); cout<<"binCenter = "<SetLineColor(4); hMassX2->SetLineColor(7); hMassX->GetXaxis()->SetRangeUser(1.2,2.4); if (iShow == -1) { hMassX->Draw(""); hMassX1->Draw("samehist"); hMassX2->Draw("samehist"); hMassXB->SetLineColor(2); hMassXB->Draw("samehist"); } if (iShow == -2){ hMassX1->SetMinimum(0); hMassX1->SetMaximum(50); hMassX1->Draw(""); } //if (iShow == -3){ //hMassX2->SetMinimum(0); //hMassX2->SetMaximum(500); //hMassX2->Draw("hist"); //hMassX2->Draw("hist"); //} if (iShow == -3) hMassX2->Draw("hist"); if (iShow == -4) hMassXB->Draw(""); if (iShow == -10) { hMassXP0->SetMaximum(400); hMassXP0->SetMinimum(0.0); hMassXP0->Draw(); } if (iShow == -11) { hMassXP1->SetMaximum(1.03); hMassXP1->SetMinimum(1.0); hMassXP1->Draw(); } if (iShow == -12) { //hMassXP2->SetMaximum(1.03); //hMassXP2->SetMinimum(1.0); hMassXP2->Draw(); } if (iShow == -13) { //hMassXP3->SetMaximum(1.03); //hMassXP3->SetMinimum(1.0); hMassXP3->Draw(); } if (iShow == -14) { //hMassXP4->SetMaximum(1.03); //hMassXP4->SetMinimum(1.0); hMassXP4->Draw(); } if (iShow == -15) { //hMassXP5->SetMaximum(1.03); //hMassXP5->SetMinimum(1.0); hMassXP5->Draw(); } if (iShow == -16) { hMassXP6->SetMaximum(8.0); hMassXP6->SetMinimum(0.0); hMassXP6->Draw(); } if (iShow == -17) { //hMassXP7->SetMaximum(1.03); //hMassXP7->SetMinimum(1.0); hMassXP7->Draw(); } if (iShow == -18) { hMassXP8->SetMaximum(0.8); hMassXP8->SetMinimum(-0.2); hMassXP8->Draw(); } if (iShow == -19) { hMassXP9->SetMaximum(1.02); hMassXP9->SetMinimum(0.97); hMassXP9->Draw(); } if (iShow == -110) { hMassXP10->SetMaximum(800.0); hMassXP10->SetMinimum(50.0); hMassXP10->Draw(); } }