Просмотр исходного кода

Added Gamma-based fitting procedure, small fixes

PeterParfenov лет назад: 4
Родитель
Сommit
56cac52c5b
3 измененных файлов с 340 добавлено и 3 удалено
  1. 4 1
      CMakeLists.txt
  2. 98 2
      FitCent.cpp
  3. 238 0
      FitCentGamma.cpp

+ 4 - 1
CMakeLists.txt

@@ -27,4 +27,7 @@ set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -O3")
 
 #---Generate executable
 add_executable(FitCent FitCent.cpp)
-target_link_libraries(FitCent ${ROOT_LIBRARIES})  
+target_link_libraries(FitCent ${ROOT_LIBRARIES})
+
+add_executable(FitCentGamma FitCentGamma.cpp)
+target_link_libraries(FitCentGamma ${ROOT_LIBRARIES})

+ 98 - 2
FitCent.cpp

@@ -43,6 +43,24 @@ Double_t ftSigmN(Double_t *x, Double_t *par)
     return func;
 }
 
+Double_t ftEtaN(Double_t *x, Double_t *par)
+{
+    // Parameters
+    Double_t sig0 = par[0];
+    Double_t n0 = par[1];
+    Double_t a1 = par[2];
+    Double_t a2 = par[3];
+    Double_t a3 = par[4];
+    // Variables
+    Double_t xx = x[0];
+    // Function
+    Double_t mean_n = ftMeanN(&x[0], &par[1]);
+    Double_t sigma  = ftSigmN(&x[0], &par[0]);
+    Double_t func = 2 * TMath::Power(1+TMath::Erf(mean_n/(TMath::Sqrt(2)*sigma)), -1);
+
+    return func;
+}
+
 Double_t ftProbIntegrand(Double_t *x, Double_t *par)
 {
     // Parameters
@@ -63,6 +81,27 @@ Double_t ftProbIntegrand(Double_t *x, Double_t *par)
     return func;
 }
 
+Double_t ftProbIntegrandEta(Double_t *x, Double_t *par)
+{
+    // Parameters
+    Double_t n = par[0];
+    Double_t sig0 = par[1];
+    Double_t n0 = par[2];
+    Double_t a1 = par[3];
+    Double_t a2 = par[4];
+    Double_t a3 = par[5];
+    // Variables
+    Double_t xx = x[0];
+    // Function
+    Double_t mean_n = ftMeanN(&x[0], &par[2]);
+    Double_t sigma = ftSigmN(&x[0], &par[1]);
+    Double_t eta   = ftEtaN(&x[0], &par[1]);
+    Double_t func = eta / (sigma * TMath::Sqrt(2 * TMath::Pi())) *
+                    TMath::Exp(-TMath::Power((n - mean_n), 2) / (2 * sigma * sigma));
+
+    return func;
+}
+
 Double_t ftProbN(Double_t *x, Double_t *par)
 {
     // Parameters
@@ -81,6 +120,24 @@ Double_t ftProbN(Double_t *x, Double_t *par)
     return func;
 }
 
+Double_t ftProbNEta(Double_t *x, Double_t *par)
+{
+    // Parameters
+    Double_t sig0 = par[0];
+    Double_t n0 = par[1];
+    Double_t a1 = par[2];
+    Double_t a2 = par[3];
+    Double_t a3 = par[4];
+    // Variables
+    Double_t n = x[0];
+    // Function
+    TF1 *f = new TF1("f", ftProbIntegrandEta, 0, 1, 6);
+    f->SetParameters(n, sig0, n0, a1, a2, a3);
+    Double_t func = f->Integral(0, 1);
+
+    return func;
+}
+
 Double_t ftCentBIntegrand(Double_t *x, Double_t *par)
 {
     // Parameters
@@ -139,6 +196,27 @@ Double_t ftProbCb(Double_t *x, Double_t *par)
     return func;
 }
 
+Double_t ftProbB(Double_t *x, Double_t *par)
+{
+    // Parameters
+    Double_t sigInel = par[0];
+    Double_t n = par[1];
+    Double_t sig0 = par[2];
+    Double_t n0 = par[3];
+    Double_t a1 = par[4];
+    Double_t a2 = par[5];
+    Double_t a3 = par[6];
+    // Variables
+    Double_t b  = x[0];
+    Double_t cb = TMath::Pi()*TMath::Power(b,2)/sigInel;
+    // Function
+    TF1 *f = new TF1("f", ftProbCb, 0, 1, 6);
+    f->SetParameters(n, sig0, n0, a1, a2, a3);
+    Double_t func = 2*TMath::Pi()*b/sigInel*f->Eval(cb);
+
+    return func;
+}
+
 int main(int argc, char **argv)
 {
     if (argc != 5)
@@ -268,7 +346,14 @@ int main(int argc, char **argv)
     for (int iter = 0; iter < NumOfIterations; iter++)
     {
         std::cout << "Fit iteration: " << iter << std::endl;
-        grMult->Fit(fitFunc, "R");
+        grMult->Fit(fitFunc, "RM");
+    }
+
+    // Fitting function w/ ftEtaN
+    TF1 *fitEtaFunc = new TF1("fitEtaFunc", ftProbNEta, vMult.at(0), vMult.at(vMult.size() - 1), 5);
+    for (int ipar = 0; ipar < 5; ipar++)
+    {
+        fitEtaFunc->SetParameter(ipar, fitFunc->GetParameter(ipar));
     }
 
     // Build Centrality vs N_{ch} function using fit parameters
@@ -278,17 +363,23 @@ int main(int argc, char **argv)
         centFunc->SetParameter(ipar, fitFunc->GetParameter(ipar));
     }
 
-    // Calculate probability function for 10% step centrality
+    // Calculate probability functions for 10% step centrality
     TF1 *probFunc[10];
+    TF1 *probFuncB[10];
     Double_t n_test;
+    const Double_t sigma_Inelastic = 777.; // [fm^2] for Pb+Pb
     for (int icent=0; icent<10; icent++)
     {
         probFunc[icent] = new TF1(Form("probFunc_cent%i",icent), ftProbCb, 0, 1, 6);
+        probFuncB[icent] = new TF1(Form("probFuncB_cent%i",icent), ftProbB, 0, 20, 7);
         n_test = centFunc->GetX(icent/10.+0.05);
         probFunc[icent]->SetParameter(0,n_test);
+        probFuncB[icent]->SetParameter(0,sigma_Inelastic);
+        probFuncB[icent]->SetParameter(1,n_test);
         for (int ipar=0; ipar<5; ipar++)
         {
             probFunc[icent]->SetParameter(ipar+1,fitFunc->GetParameter(ipar));
+            probFuncB[icent]->SetParameter(ipar+2,fitFunc->GetParameter(ipar));
         }
     }
 
@@ -300,11 +391,16 @@ int main(int argc, char **argv)
 
     grMult->Write();
     fitFunc->Write();
+    fitEtaFunc->Write();
     centFunc->Write();
     for (int icent=0; icent<10; icent++)
     {
         probFunc[icent]->Write();
     }
+    for (int icent=0; icent<10; icent++)
+    {
+        probFuncB[icent]->Write();
+    }
 
     fo->Close();
 

+ 238 - 0
FitCentGamma.cpp

@@ -0,0 +1,238 @@
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include <string>
+#include <TROOT.h>
+#include <TMath.h>
+#include <TF1.h>
+#include <TH1F.h>
+#include <TGraphErrors.h>
+#include <TFile.h>
+#include <TString.h>
+#include <TStopwatch.h>
+#include <Math/PdfFuncMathCore.h>
+#include <Math/IntegratorOptions.h>
+
+Double_t ftKcb(Double_t *x, Double_t *par)
+{
+    // Parameters
+    Double_t theta = par[0];
+    Double_t k_max = par[1];
+    Double_t a1 = par[2];
+    Double_t a2 = par[3];
+    Double_t a3 = par[4];
+    Double_t a4 = par[5];
+    // Vartiables
+    Double_t xx = x[0];
+    // Function
+    Double_t func = k_max/theta * TMath::Exp(a1 * xx + a2 * xx * xx + a3 * xx * xx * xx + a4 * xx * xx * xx * xx);
+
+    return func;
+}
+
+Double_t ftPnIntegrand(Double_t *x, Double_t *par)
+{
+    // Parameters
+    Double_t n = par[0];
+    Double_t theta = par[1];
+    Double_t k_max = par[2];
+    Double_t a1 = par[3];
+    Double_t a2 = par[4];
+    Double_t a3 = par[5];
+    Double_t a4 = par[6];
+    // Vartiables
+    Double_t xx = x[0];
+    // Function
+    Double_t k = ftKcb(&x[0], &par[1]);
+    Double_t func = ROOT::Math::gamma_pdf(n,k,theta);
+
+    return func;
+}
+
+Double_t ftPn(Double_t *x, Double_t *par)
+{
+    // Parameters
+    Double_t theta = par[0];
+    Double_t k_max = par[1];
+    Double_t a1 = par[2];
+    Double_t a2 = par[3];
+    Double_t a3 = par[4];
+    Double_t a4 = par[5];
+    // Variables
+    Double_t n  = x[0];
+    // Function
+    TF1 *f = new TF1("f",ftPnIntegrand,0,1,7);
+    f->SetParameters(n,theta,k_max,a1,a2,a3,a4);
+    Double_t func = f->Integral(0,1);
+
+    return func;
+}
+
+Double_t ftCent(Double_t *x, Double_t *par)
+{
+    // Parameters
+    Double_t theta = par[0];
+    Double_t k_max = par[1];
+    Double_t a1 = par[2];
+    Double_t a2 = par[3];
+    Double_t a3 = par[4];
+    Double_t a4 = par[5];
+    // Variables
+    Double_t n  = x[0];
+    // Function
+    TF1 *f = new TF1("f",ftPn,0,500);
+    f->SetParameters(theta,k_max,a1,a2,a3,a4);
+    Double_t func = f->Integral(n,TMath::Infinity());
+
+    return func;
+}
+
+int main(int argc, char **argv)
+{
+    if (argc != 5)
+    {
+        std::cerr << "Error: wrong number of arguments:" << std::endl;
+        std::cerr << "./FitCentGamma -dat/-root inputfile -o output" << std::endl;
+        return 10;
+    }
+
+    TStopwatch timer;
+    timer.Start();
+
+    ROOT::Math::IntegratorOneDimOptions::SetDefaultRelTolerance(1.E-3); 
+
+    // Number of iterations during fitting procedure
+    int NumOfIterations = 5;
+
+    // Parsing arguments
+    TString inputfile, outputfile;
+    Bool_t isRoot = false;
+    for (int i = 1; i < argc; i++)
+    {
+        if (std::string(argv[i]) != "-dat" &&
+            std::string(argv[i]) != "-root" &&
+            std::string(argv[i]) != "-o")
+        {
+            std::cerr << "Error: Unknown option: " << argv[i] << std::endl;
+            return 11;
+        }
+        else
+        {
+            if ((std::string(argv[i]) == "-dat" || std::string(argv[i]) == "-dat") && i == argc - 1)
+            {
+                std::cerr << "Error: input file was not specified!" << std::endl;
+                return 12;
+            }
+            if (std::string(argv[i]) == "-root" && i != argc - 1)
+            {
+                inputfile = argv[++i];
+                isRoot = true;
+            }
+            if (std::string(argv[i]) == "-dat" && i != argc - 1)
+            {
+                inputfile = argv[++i];
+                isRoot = false;
+            }
+            if (std::string(argv[i]) == "-o" && i == argc - 1)
+            {
+                std::cerr << "Error: output file was not specified!" << std::endl;
+                return 13;
+            }
+            if (std::string(argv[i]) == "-o" && i != argc - 1)
+            {
+                outputfile = argv[++i];
+            }
+        }
+    }
+
+    std::cout << "Input file:  " << inputfile << std::endl;
+    std::cout << "Output file: " << outputfile << std::endl;
+    std::vector<Double_t> vMult, vProb, vError, vErrorX;
+    Double_t mult, prob, error;
+    int counter = 0;
+    if (!isRoot)
+    {
+        std::ifstream infile(inputfile.Data());
+        if (!infile.is_open())
+        {
+            std::cerr << "Error: can not open the file!" << std::endl;
+            return 20;
+        }
+        std::string str;
+        const Double_t bin_width = 0.44;
+        while (std::getline(infile, str))
+        {
+            std::istringstream iss(str);
+            if (!(iss >> mult >> prob))
+            {
+                std::cerr << "Error: can not read line " << counter << ": " << str << std::endl;
+                return 30;
+            }
+
+            error = prob / (1e7 * 0.44);
+
+            vMult.push_back(mult);
+            vProb.push_back(prob);
+            vError.push_back(error);
+            vErrorX.push_back(0.);
+
+            counter++;
+        }
+        std::cout << counter << " lines were processed." << std::endl;
+    }
+    if (isRoot)
+    {
+        TFile *fi = new TFile(inputfile.Data(), "read");
+        TH1F *hist = (TH1F *)fi->Get("hRefMultSTAR");
+        if (!hist)
+        {
+            std::cerr << "Could not read histogram from root-file!" << std::endl;
+            return 31;
+        }
+
+        for (int i = 0; i < hist->GetNbinsX(); i++)
+        {
+            if (hist->GetBinContent(i + 1) == 0)
+                continue;
+            vMult.push_back(hist->GetBinCenter(i + 1));
+            vProb.push_back(hist->GetBinContent(i + 1) / hist->Integral());
+            vError.push_back(hist->GetBinError(i + 1) / hist->Integral());
+            vErrorX.push_back(0.);
+            counter++;
+        }
+        std::cout << counter << " bins were processed." << std::endl;
+        delete hist;
+        fi->Close();
+    }
+
+    // Getting input data into TGraphErrors
+    TGraphErrors *grMult = new TGraphErrors(vMult.size(), &vMult[0], &vProb[0], &vErrorX[0], &vError[0]);
+    // TGraph *grMult = new TGraph(vMult.size(), &vMult[0], &vProb[0]);
+    grMult->SetName("grMult");
+    grMult->GetYaxis()->SetTitle("1/N dN/dN_{ch}");
+    grMult->GetXaxis()->SetTitle("N_{ch}");
+
+    // Fitting procedure
+    TF1 *fitFunc = new TF1("fitFunc", ftPn, vMult.at(0), vMult.at(vMult.size() - 1), 6);
+    fitFunc->SetParameters(0.3,183.,-4.,3.,-7.,2.);
+
+    for (int iter = 0; iter < NumOfIterations; iter++)
+    {
+        grMult->Fit(fitFunc,"RM");
+    }
+
+    std::cout << std::endl;
+    std::cout << "Writing output." << std::endl;
+
+    TFile *fo = new TFile(outputfile.Data(), "recreate");
+    fo->cd();
+
+    grMult->Write();
+    fitFunc->Write();
+
+    fo->Close();
+
+    timer.Stop();
+    timer.Print();
+    return 0;
+}