-
-//______________________________________________
-void AliDielectronSignalFunc::ProcessFitIKF(TObjArray * const arrhist) {
- //
- // Fit the +- invariant mass distribution only
- //
- //
-
- const Double_t bigNumber = 100000.;
- Double_t chi2ndfm1 = bigNumber;
- Double_t ratiom1 = bigNumber;
- Double_t chi2ndf = bigNumber;
- Int_t nDP =0;
-
- Int_t maxPolDeg = 8;
-
- fHistDataPM = (TH1F*)(arrhist->At(1))->Clone("histPM"); // +- SE
- if(fRebin>1) fHistDataPM->Rebin(fRebin);
-
- fgHistSimPM = (TH1F*)(arrhist->At(3))->Clone("histPMsim"); // +- mc shape
- if (!fgHistSimPM) {
- AliFatal("No mc peak shape found at idx 3.");
- return;
- }
- if(fRebin>1) fgHistSimPM->Rebin(fRebin);
-
- // try out the polynomial degrees
- for (Int_t iPD=0; iPD<=maxPolDeg; iPD++) {
- TH1F *hf1 = (TH1F *) fHistDataPM->Clone(Form("hf1_PD%d",iPD));
-
- FitOneMinv(hf1, fgHistSimPM, iPD);
- if (fChi2Dof > 0) chi2ndf = fChi2Dof;
- AliInfo(Form("nDP: %d, iPD: %d, chi2ndf: %f", nDP, iPD, chi2ndf));
-
- ratiom1 = TMath::Abs(fChi2Dof - 1);
- if (chi2ndfm1 > ratiom1) { // search for the closest to 1.
- chi2ndfm1 = ratiom1;
- nDP = iPD;
- }
- }
-
-
- // fit again with the best polynomial degree
- TH1F *h2 = (TH1F *) fHistDataPM->Clone(Form("h2_PD%d",nDP));
-
- FitOneMinv(h2, fgHistSimPM, nDP);
- AliInfo(Form("Best Fit: PD %d, chi^2/ndf %.3f, S/B %.3f",nDP,fChi2Dof,fValues(3)));
-
-}
-//______________________________________________
-void AliDielectronSignalFunc::FitOneMinv(TH1F *hMinv, TH1F *hSim, Int_t pod) {
- //
- // main function to fit an inv mass spectrum
- //
-
- TObjArray *arrResults = new TObjArray;
- arrResults->SetOwner();
- arrResults->AddAt(hMinv,0);
-
- // Degree of polynomial
- fPolDeg = pod;
-
- // inclusion and exclusion areas (values)
- const Double_t kJPMass = 3.096916;
- // inclusion and exclusion areas (bin numbers)
- const Int_t binIntLo = hMinv->FindBin(fIntMin);
- const Int_t binIntHi = hMinv->FindBin(fIntMax);
- // for error calculation
- Double_t intAreaEdgeLo = hMinv->GetBinLowEdge(binIntLo);
- Double_t intAreaEdgeHi = hMinv->GetBinLowEdge(binIntHi)+hMinv->GetBinWidth(binIntHi);
- Double_t norm = (binIntHi-binIntLo)/(fIntMax - fIntMin);
-
- TH1F *hBfit = (TH1F *) hMinv->Clone(); // for bg only fit (excluding peak region)
- TH1F *hSigF = (TH1F *) hMinv->Clone(); // signal with subtracted bg
- TH1F *hBgrd = (TH1F *) hMinv->Clone(); // bg histogram
-
- hBfit->Reset();
- hSigF->Reset();
- hBgrd->Reset();
-
- // extract start parameter for the MC signal fit
- Double_t bgvalAv = (hMinv->Integral(1,hMinv->GetNbinsX()+1) - hMinv->Integral(binIntLo,binIntHi)) / (hMinv->GetNbinsX()+1 - (binIntHi-binIntLo));
- Double_t pkval = hMinv->GetBinContent(hMinv->FindBin(kJPMass)) - bgvalAv;
- Double_t heightMC = hSim->GetBinContent(hSim->FindBin(kJPMass));
- Double_t peakScale = (heightMC > 0. ? pkval/heightMC : 0.0);
-
- Int_t nBgnd = 2 + fPolDeg; // degree + c1st oefficient + higher coefficients
- Int_t nMinv = nBgnd + 1; // bgrd + peakscale
-
- // Create the spectra without peak region
- for (Int_t iBin = 0; iBin <= hMinv->GetNbinsX(); iBin++) {
- if ((iBin < binIntLo) || (iBin > binIntHi)) {
- hBfit->SetBinContent(iBin,hMinv->GetBinContent(iBin));
- hBfit->SetBinError(iBin,hMinv->GetBinError(iBin));
- }
- }
-
-
- // =======
- // 1.
- // =======
- // Do the fit to the background spectrum
- TF1 *fBo = new TF1("bgrd_fit",BgndFun,fFitMin,fFitMax,nBgnd);
- for (Int_t iPar=0; iPar<nBgnd; iPar++) {
- if (iPar == 0) fBo->FixParameter(0, fPolDeg);
- if (iPar == 1) fBo->SetParameter(iPar, bgvalAv);
- if (iPar >= 2) fBo->SetParameter(iPar, 0.);
- }
- hBfit->Fit(fBo,"0qR");
- //hBfit->SetNameTitle("bgrd_fit");
- arrResults->AddAt(fBo,1);
-
-
- // =======
- // 2.
- // =======
- // Fit the whole spectrum with peak and background
- TF1 *fSB = new TF1("bgrd_peak_fit",MinvFun,fFitMin,fFitMax,nMinv);
- fSB->FixParameter(0, fPolDeg);
- fSB->SetParameter(1, peakScale);
- // copy the polynomial parameters
- for (Int_t iPar=0; iPar<=fPolDeg; iPar++)
- fSB->SetParameter(2+iPar, fBo->GetParameter(iPar+1));
- hMinv->Fit(fSB,"0qR");
- arrResults->AddAt(fSB,2);
-
-
- // =======
- // 3.
- // =======
- // Create the background function
- TF1 *fB = new TF1("bgrdOnly_fkt",BgndFun,fFitMin,fFitMax,nBgnd);
- fB->FixParameter(0,fPolDeg);
- for (Int_t iDeg=0; iDeg<=fPolDeg; iDeg++) {
- fB->SetParameter(1+iDeg,fSB->GetParameter(2+iDeg));
- fB->SetParError(1+iDeg,fSB->GetParError(2+iDeg));
- }
- // create background histogram from background function
- hBgrd->Eval(fB);
- hBgrd->Fit(fB,"0qR");
- // calculate the integral and integral error from fit function
- Double_t intc = fB->Integral(intAreaEdgeLo, intAreaEdgeHi) * norm;
- Double_t inte = fB->IntegralError(intAreaEdgeLo, intAreaEdgeHi) * norm;
- arrResults->AddAt(fB,3);
-
- // Fill the background spectrum erros. Use the error from the fit function for the background fB
- for (Int_t iBin = 0; iBin <= hBgrd->GetNbinsX(); iBin++) {
- Double_t x = hBgrd->GetBinCenter(iBin);
- if ((x >= fFitMin) && (x <= fFitMax)) {
- Double_t binte = inte / TMath::Sqrt((binIntHi-binIntLo)+1);
- hBgrd->SetBinError(iBin,binte);
- }
- }
- arrResults->AddAt(hBgrd,4);
-
- // =======
- // 4.
- // =======
- // Subtract the background
- hSigF->Add(hMinv,hBgrd,1.0,-1.0);
- for (Int_t iBin = 0; iBin <= hSigF->GetNbinsX(); iBin++) {
- Double_t x = hSigF->GetBinCenter(iBin);
- if ((x < fFitMin) || (x > fFitMax)) {
- hSigF->SetBinContent(iBin,0.0);
- hSigF->SetBinError(iBin,0.0);
- }
- }
- hSigF->SetNameTitle("peak_only","");
- arrResults->AddAt(hSigF,5);
-
- // =======
- // 5.
- // =======
- // Fit the background-subtracted spectrum
- TF1 *fS = new TF1("peakOnly_fit",PeakFunCB,fFitMin,fFitMax,5);
- fS->SetParameters(-.05,1,kJPMass,.003,700);
- fS->SetParNames("alpha","n","meanx","sigma","N");
- hSigF->Fit(fS,"0qR");
- arrResults->AddAt(fS,6);
-
-
- // connect data members
- fFuncSignal = (TF1*) arrResults->At(6)->Clone();
- fFuncBackground = (TF1*) arrResults->At(3)->Clone();
- fFuncSigBack = (TF1*) arrResults->At(2)->Clone();
- fHistSignal = (TH1F*)arrResults->At(5)->Clone();
- fHistBackground = (TH1F*)arrResults->At(4)->Clone();
-
- // fit results
- Double_t chi2 = fSB->GetChisquare();
- Int_t ndf = fSB->GetNDF();
- fChi2Dof = chi2/ndf;
-
- // signal + signal error
- fValues(0) = hSigF->IntegralAndError(binIntLo, binIntHi, fErrors(0));
- fValues(1) = intc; // background
- fErrors(1) = inte; // background error
- // S/B (2) and significance (3)
- SetSignificanceAndSOB();
- fValues(4) = fS->GetParameter(2); // mass
- fErrors(4) = fS->GetParError(2); // mass error
- fValues(5) = fS->GetParameter(3); // mass wdth
- fErrors(5) = fS->GetParError(3); // mass wdth error
-
-
- delete arrResults;
-
-}
-//______________________________________________________________________________
-Double_t AliDielectronSignalFunc::BgndFun(const Double_t *x, const Double_t *par) {
- // parameters
- // [0]: degree of polynomial
- // [1]: constant polynomial coefficient
- // [2]..: higher polynomial coefficients
-
- Int_t deg = ((Int_t) par[0]);
-
- Double_t f = 0.0;
- Double_t yy = 1.0;
- Double_t xx = x[0];
-
- for (Int_t i = 0; i <= deg; i++) {
- f += par[i+1] * yy;
- yy *= xx;
- }
-
-
- return f;
-}