W opisanym tu ćwiczeniu zajmę się najprostszym fitowaniem - fitowaniem 1D funkcji, której przebieg analityczny jest znany. Najpierw musimy w naszym makrze stworzyć jakąś funkcję korelacyjną żeby było co dopasowywać.
Hal::Femto1DCF* cf = new Hal::Femto1DCF("func", 100, 0, 1, Hal::Femto::EKinematics::kPRF);
// fill dummy data
for (int i = 1; i <= 100; i++) {
double kstar = cf->GetNum()->GetBinCenter(i);
cf->GetNum()->SetBinContent(i, 1.0 + TMath::Exp(-400. * kstar * kstar));
cf->GetNum()->SetBinError(i, 0.01);
cf->GetDen()->SetBinContent(i, 10);
cf->GetDen()->SetBinError(i, 0.01);
}
cf->GetNum()->SetMarkerStyle(kFullCircle);
Warto zwrócić uwagę na ostatnią linię - manipulując opcjami rysowania histogramu licznika (styl, kolor markerów itd.) manipulujemy jednocześnie style rysowanej funkcji korelacyjnej.
Kiedy na szybko stworzyliśmy sobie funkcję korelacyjną odpowiadającą z grubsza R=2fm. Teraz budujemy sobie klasę fitującą, Tutaj ustawiamy zakres parametrów - działa to podobnie jak w TF1 mamy FixParameter czy SetParLimits, ustawiamy sobie też zakres dopasowania.
auto fit = new Hal::CorrFit1DCF_Gauss(); fit->SetParLimits(fit->RadiusID(), 1, 5); fit->SetRLimits(1, 5); // alias do funkcji wyzej fit->SetParLimits(fit->LambdaID(), 0.5, 1); fit->SetLambdaLimits(0.5, 1); // alias do funkcji wyżej fit->FixParameter(fit->NormID(), 0.1); fit->SetRange(0, 0.75);//zakres dopasowania
Następnie zaczynają się istotne różnice. Przy fitowaniu funkcji korelacyjnych jest kilka istotnych problemów.
Pierwszy - nie zawsze chcemy fitować "po zwykłym" chi2, dlatego corrfit HAL-owy posiada 3 rodzaje funkcji minimalizowanych (SetMinimzedFunc): kChi minimalizuje funkcję chi2 między korelacyjną pomnożoną przez mianownik fitowaną do licznika czyli coś takiego:
kChi2 to standardowa minimalizacja chi2 między funkcją korelacyjną teoretyczną, a histogramem z eksperymentu, trzecia opcja kLog - to minimalizacja zalecana przy małej liczbie wejść w danych.
Drugi problem to, że ROOT ma problemy ze znalezieniem funkcji, to powoduje, że często fity są po prostu złe. Tutaj zastosowałem pewne obejście - zamiast używać zwykłego ROOTowego minimizera (np. EMinAlgo::kMinuitMigrad który używa tylko minuita) używam dwustopniowego minimizera - kHalScanMigrad). Jak on działa? używa pliku XML który wygląda np. tak:
<minimizer> <param name="N" min="0" max="1" step="0.01"></param> <param name="R" min="1" max="10" step="0.01"></param> <param name="#lambda" min="0" max="1" step="0.01"></param> </minimizer>
Ważne: plik musi określać liczbę kroków wszystkich parametrów, nawet jeśli są ustawione na sztywno, a zakresy powinny być co najmniej tak duże jak zakresy które chcemy ustawić przy dopasowywaniu!
HAL teraz porównuje sobie plik XML z ustawieniami (SetParLimits) i testuje wszystkie dozwolone kombinacje parametrów z zadanymi krokami. Gdy znajdzie minimum dopasowania przełącza się na ROOTowy minimizer. Teraz parametrem startowym jest punkt gdzie znaleziono minimum, a dozwolone zakresy dopasowania to +/- jeden punkt (chyba, że parametr jest zafiksowany). Im więcej punktów, tym oczywiście większa szansa, że fit będzie poprawny (bo trafimy bliżej globalnego minimum i ROOT-owy minimizer będzie miał łatwiej) - ale czas fitu szybko rośnie wraz z liczbą punktów i liczbą wolnych parametrów.
auto conf = Hal::MinimizerStepConf();
conf.LoadFromXML("conf.xml");
fit->SetMinimizerConf("conf.xml");
fit->SetLineColor(kRed + 2);
fit->SetMinimizedFunc(Hal::CorrFit::EMinFunc::kChi2);
fit->SetMinimizer(Hal::CorrFit::EMinAlgo::kHalScanMigrad); // kMInimizerScan = hardorowy minimizer
// TCanvas* c = new TCanvas();
cf->Fit(fit);
Możemy też funkcję narysować wraz z dopasowaniem
TCanvas* chi2 = new TCanvas();
chi2->Divide(1, 2);
chi2->cd(1);
cf->Draw("grid+{y=0,2}+pad");
fit->Draw("legend+chi2s+norm");
Zauważyć tu można następujące opcje rysowania eksperymentalnej funkcji korelacyjnej:
- pad - rysuj na aktualnym padzie (domyślnie HAL tworzy własny)
- grid - rysuj siatkę na padzie
- {y=0,2} - zakres OY 0-2
W przypadku funkcji dopasowanej mamy opcje:
- legend - rysuj legendę dopasowania
- chi2s - rysuj chi2 z dopasowania na legendzie
- norm - po narysowaniu przeskaluj rysowaną funkcję korelacyjną, tak żeby teoretyczna saturowała się na 1
Na końcu sobie rysujemy mapę chi2 na której rysujemy chi2 w funkcji promienia (RadiusID) i lambdy (LambdaID)
chi2->cd(2);
ChiSqMap2D* chimap = fit->GetChiSquareMap(fit->RadiusID(), 200, 1, 10, fit->LambdaID(), 200, 0.5, 1.5);
chimap->Draw("min");
Opcjonalnie możemy sobie uruchomić GUI gdzie sobie możemy wyklikać dopasowanie.
Hal::CorrFitGUI* gui = new Hal::CorrFitGUI(fit);